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ABSTRACT 


A computer program was developed and written in Fortran V which predicts the 
transient and steady state performance and heat transfer characteristics of a pulsing 
G0 2 /GH 2 rocket engine. This program predicts the dynamic flow and ignition charac- 
teristics which, when combined in a quasi-steady state manner with the combustion and 
mixing analysis program, will provide the thrust and specific impulse of the engine as 
a function of time. The program also predicts the transient and steady state heat trans- 
fer characteristics of the engine using various cooling concepts. The computer program, 
test case , and documentation are presented in this document. The program is applicable 
to any system capable of utilizing the Fortran IV or Fortran V language. 
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FOREWORD 


Contract NAS 9-11740, "Characterization of the Space Shuttle Reaction Control 
System Engine" was performed by The Marquardt Company at Van Nuys , California, 
and by the principal Subcontractor, The General Applied Science Laboratory (GASL) 
The period of performance covered by this report is from May 15, 1971, to April 15, 
1972. The Program Manager was R. C. Stechman. Principal Investigators included 
M. S. Wilson, R. B. Edelman, and O. Fortune. The NASA Project Manager was 
D. Hyatt, NASA-MSC. 
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CHARACTERIZATION OF THE SPACE SHUTTLE 
REACTION CONTROL SYSTEM ENGINE 


By: M. S. Wilson, R. C. Stechman, R. B. Edelman, 
O. F. Fortune, and C. Economos 

The Marquardt Company 
Van Nuys , California 


SUMMARY 


Under Contract NAS 9-11740, a computer model has been developed to 
characterize the combustion, ignition, and dynamic performance of G0 2 /GH 2 rocket 
engines. A dynamics module, incorporating valve opening characteristics and ignition 
system performance, is used to predict pressures, temperatures, and mixture ratios 
. in the various engine cavities during pulse operation. Instantaneous performance may 
be calculated at intervals utilizing a combustion model to determine the degree of 
combustion gas mixing and a performance model to calculate friction losses, diver- 
gence losses, etc. The combustion model may be used independently to predict the 
chemical and gas dynamic behavior of the reacting O 2 /H 2 gas mixture considering 
either finite or infinite chemical reaction rates and including the effects of mass 
addition from film, transpiration, or dump cooling. The characterization program 
also includes a heat transfer module designed to predict operating temperatures for 
G0 2 /GH 2 rocket engines using radiation, heatsink, film, regenerative, or dump 
(liner) cooling schemes. 

Analytical results from the characterization models have been found to com- 
pare favorably with test results obtained with G0 2 /GH 2 engines. Comparisons were 
made with data on dynamic response, specific impulse, and structural temperatures. 
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INTRODUCTION 

Development of analytical techniques and computer models which can characterize 
various parts of the space shuttle is required in order that the overall system can be 
analyzed without extensive testing of each proposed design. The models developed can 
be used to determine system design criteria, influence coefficients, and ultimately to 
characterize the final system design. An important aspect of the shuttle system is a 
model of the reaction control system engine. Included in this model is a requirement 
for a more thorough understanding of the combustion and ignition process and their coupling 
with the dynamic and heat transfer characteristics of the engine. The successful incor- 
poration of the combustion and mixing model with other parametric models provides a 
capability to determine the reaction control system performance for a given point design 
as well as a fixed set of environmental conditions . 

This report presents the results of a 10-month study for the National Aeronautics 
and Space Administration/Manned Spacecraft Center to develop an analytical model and 
computer program for the prediction of the transient and steady state performance charac- 
teristics of reaction control rocket engines which use gaseous oxygen and gaseous hydrogen. 
Major emphasis has been placed on the development of computer codes to describe the 
mixing and reaction processes in the combustion chamber. Less emphasis has been placed 
on the losses and processes in the exit nozzle since these processes have been previously 
investigated in great detail (Reference 1). The results generated by the program when com- 
pared to the results of the G 0 2 /GH 2 studies now being conducted will provide the necessary 
data for meaningful correlation. 

This report is divided into two discrete parts. The first part describes the metho- 
dology used to provide the input for the computer programming. The six parts of the 
model, dynamics, injection, ignition, combustion, performance, and heat transfer are 
described and the equations are presented in detail. The results of a comparison of a 
single test case with the data obtained from a contractor's report (Reference 2) are des- 
cribed. 

The second portion of the report contains a detailed description of the computer 
program, a listing of the program, a description of input requirements, and flow charts 
which will enable the reader to use the program. The description and results of the test 
case are also provided. 
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ANALYSIS 


Model Description 

The performance of reaction control engines is measured by the instantaneous 
thrust that the engine generates and the amount of propellant used in producing the thrust. 
The performance of G0 2 /GH2 engines is basically a function of the parameters listed in 
Table I. These parameters can, when described adequately, predict the specific impulse 
and heat transfer characteristics of an engine in both a pulsing cycle and at steady state. 

As shown, the parameters which affect specific impulse are more numerous than those that 
affect thrust and mixture ratio. Thrust variations are only a function of system changes 
and can be varied and thus are not a true evaluator of efficiency of the engine compared to 
specific impulse. Mixture ratio variations are also caused principally by system variation 
and time (gas dynamics). Therefore, in the performance analysis of an engine, specific 
impulse efficiency is the true measure of performance while thrust and mixture ratio 
characterize the engine due to system changes and time. Heat transfer characteristics, 
in terms of temperature variation, will influence performance in a minor way. 

Based on the above premise, a fully integrated pulsing performance model, based 
on variations of the system shown in Figure 1, was developed for gaseous oxygen/gaseous 
hydrogen rocket engines. The model predicts instantaneous thrust and integrated impulse . 
and Isp for specified pulsing conditions taking account for injector mixing, turbulent diffu- 
sion, finite chemical kinetics, divergence losses, friction losses, and heat loss. 

The basic driver for the program is an engine dynamics model which calculates 
pressure, mixture ratio, and temperature history for the engine main combustor, pilot 
combustor, manifolds, etc. The dynamics subprogram incorporates an integral spark 
ignitor model which tests for sparking, ignition, and flame quenching. Valve modeling 
allows for sequencing and finite opening and closing times. 

A propellant injection model provides a starting flow of specified mixing efficiency 
to a combustion model developed by General Applied Science Laboratory. The injection 
and combustion models are called a number of times during the chamber pressurization 
process. The combustion model predicts propellant mixing and chemistry from the injector 
face down to the throat. At the user's option, finite chemical kinetics are considered and 
the effects of film or transpiration cooling mass addition or regenerative cooling are 
accounted for. 

The performance model uses throat enthalpy profiles to predict quasi-steady state 
performance. The effects of nozzle divergence, boundary layer losses, and mass addition 
in the bell are included to obtain estimates of actual performance. Before ignition has 
occurred, a cold flow performance model is used to estimate thrust and Isp. The final 
performance calculations return an updated thrust coefficient value to the dynamics model 
where accumulated impulse and Isp are determined at each iteration. 
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Besides predicting performance, the computer program includes a heat transfer 
model designed to predict both steady state and transient engine temperatures. The heat 
transfer program creates a thermal model from a minimum of input information. Solution 
is obtained by finite difference techniques. Cooling options allow the user to specify single 
and multislot film cooling, regenerative cooling, liner cooling, or combinations of the 
three. 

Each of the main submodels of the engine characterization program may be run by 
itself by providing additional input. The dynamics model can be used alone to study pulsing 
operation. The combustion model can be . utilized alone to describe mixing and combustion 
of any oxygen-hydrogen ducted or free jet flow. And the heat transfer model may stand by 
itself to study engine thermal behavior. 


Dynamics . - The dynamics model is designed to predict pressures, temperatures, 
and mixture ratios in the cavities and volumes of a GOg/GHg rocket engine. A spark ignition 
model is an integral part of the analysis, and instantaneous and integrated performance 
parameters are part of the output. 

Figure 2 shows a dynamics model of an Aerojet GO 2 /GH 2 rocket engine described in 
Reference 2. The model consists of large volume oxygen and hydrogen supply tanks, in- 
jector manifolds, a pilot combustor, the. main combustor, and the space sink. The 
simulation of the valves feeding the pilot and the main combustor requires input opening 
and closing times and valve flow area is assumed to vary linearly with time. 

The flow' passages connecting mass accumulation volumes may be of two types. If 
the passage is identified as an orifice, entrance effects are assumed to dominate friction 
losses, and a discharge coefficient is specified. If friction is important, the passage is 
identified as a duct, and duct length is specified. 

Instanteneous thrust is calculated at each 'iteration based upon chamber pressure 
and input thrust coefficient. If the combustion and performance modules have been called, a 
computed thrust coefficient is passed to the dynamics model for performance computations. 
Integrated impulse and specific impulse are also calculated. 

Injection. - The injection model was written to provide starting profiles of gas 
velocity, mixture ratio, and temperature used for initialization of the combustion module. 

The model assumes a flow field consisting of concentric annuli of uniform flow alternating 
between plugs which are oxygen rich and those which are hydrogen rich relative to the over- 
all mixture ratioT "The^ number of distinct annuli is a function of the number of injection 
elements . ' 

The velocity, temperature, and flow area for the hydrogen rich and oxygen rich 
regions are determined to satisfy continuity, momentum, and energy, and to simulate the 
mixing produced by the injection elements. The degree of injector mixing is specified by 
an input mixing factor identical to the one usually used to characterize experimental injector 
performance . 
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The mass flow, overall mixture ratio, and total pressure required as input by the 
injection model are provided by the dynamics module when transient combustion and per- 
formance calculations are made. Otherwise, these parameters are user specified. 

Ignition. - The ignition model is an integral part of the dynamics subroutine. Its 
purpose is to test for the occurrence of ignition each time a voltage surge is applied to the 
spark plug. The dynamics model provides the instantaneous mixture ratio and pressure in 
the vicinity of the plug. When ignition is indicated, the combustor or pilot pressure and 
temperature are assumed to rise instantly to conditions of chemical equilibrium. 

In engine models which incorporate a pilot combustor, the spark plug is assumed to 
be in the pilot. Otherwise it is assumed to be in the combustor'. The ignition model checks 
first for the occurrence of a spark based upon the spark gap, the local pressure, and the 
applied potential. If a spark occurs, an ignition test is made which depends upon spark 
energy and local mixture ratio. If ignition occurs and the model includes a pilot, a test is 
made for flame quenching in the combustor. The ignition model is based upon theoretical 
considerations and upon a number of recent C^/Hg ignition test results. 

Combustion. - The high performance demands placed upon the proposed space 
shuttle systems has required that more sophisticated prediction techniques be developed for 
the detailed analysis of the thrust chamber combustion and mixing processes. Most per- 
formance evaluations are based upon one -dimensional models. These models may provide 
adequate definition of the potential performance of a particular motor, but can, at best, 
account for non-ideal behavior by using correction factors based upon existing engine data. 
This method may be sufficient for engines of similar, well investigated design, but relatively 
new thrust chamber designs are not amenable to this type of treatment. 

The processes occurring within the thrust chamber are controlled by mechanisms 
involving fluid mechanical and chemical kinetic phenomena. Turbulent mixing, reaction 
kinetics and mass transfer along the chamber walls contribute to the overall performance 
of the thrust chamber. To date little has been done to quantify the coupling of these processes 
within rocket thrust chambers. 

Accordingly, the combustion module was developed using a rather detailed mathe- 
matical model designed to describe the turbulent reacting thrust chamber flow field including 
the effects of film and transpiration cooling. 

Performance. - The principal purpose of the performance model is to determine 
the thrust and specific impulse of the engine at a specified point in time (the time at which a 
com bus ti"dn“ calculation is~ m ade ) . The enthalpy and mixture ratio at each grid point of the 
combustion program is used as the starting point for the kinetic recombination in the nozzle 
using a modified Bray criteria. The Isp at each streamline is then mass averaged over all 
streamlines to find the overall Isp. Boundary layer and divergence losses are deducted 
along with the Isp loss due to mass addition in the exhaust nozzle. 
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Heat Transfer. - The heat transfer model is designed to predict engine structure 
and coolant temperature distribution for G0 2 /GH 2 engines utilizing a number of cooling 
options. The model exchanges no information with other parts of the characterization pro- 
gram and, hence, can be run independently. 

At the user's option, the heat transfer subprogram will model the following cooling 
techniques : 

(1) Conduction cooling J 

(2) Heat sink cooling \ Passive methods 

(3) Radiation cooling 1 

(4) Film cooling 1 

(5) Regenerative cooling j Active methods 

(6) Combustor liner cooling J 

Film cooling, with either single or multiple injection stations, may be used with any of the 
other methods. 

In general, the user is required only to specify the combustor shape and combustor 
material characteristics plus coolant mass flow and delivery geometry. It is assumed that 
hydrogen gas is used for all active cooling. 

The heat transfer subroutine has provision for including an injector thermal model. 
When this option is utilized, the heat transfer characteristics of the injector must be 
specified. The subroutine computes transient or steady state injector temperature using 
the input characteristics and accounting for conduction and radiation from the combustor. 

The heat transfer program creates from the input a thermal network which is solved 
by a finite difference method. All thermal admittances are automatically computed. 


Analysis Methods 

Dynamics Analysis . - The dynamics model predicts transient and steady state 
pressures, mixture ratios, and temperatures in the various gas accumulation volumes 
of a GH 2 /G0 2 rocket engine. The gas accumulation volumes are the manifolds, regenera- 
tive cooling passages, the combustion chamber, the pilot combustor, etc. , which are 
connected by a system of flow resistances (orifices and ducts). A typical system is shown 
in Figure 2. 

The dynamics model is solved using a finite difference technique. At a given instant 
in time, the mass flows between the volumes are calculated. Then for a small increment 
of time, a new system pressure distribution is computed along with new temperatures and 
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mixture ratios. An integral spark ignition model tests for ignition at sparking times. 

Point and integrated performance parameters are calculated using a thrust coefficient 
passed from the performance model. 

Mass flow calculations: Two types of flow restrictions are provided in the dynamics 
model. When friction is unimportant, an orifice model is used requiring specification of the 
discharge coefficient. A duct model is also available for flow restrictions which are charac- 
terized by friction losses that are large relative to entrance losses. For either case, a 
flow admittance is calculated conforming to the following definition 


* = A d( P 0! - Po 2 ) 


For orifices, the admittance is given by 
AC ,P i 


= ( p °i - p ^r 

' (y) for choked 


flow 




The choking pressure ratio is given by 




choking 




y~i 


y tl 


(i) 


(2) 


(3) 


for unchoked flow 


(4) 


For ducts, where pressure loss due to friction is important, mass flow depends on the 
friction factor parameter, 4£b/^, as well as the pressures. The choking pressure 
ratio is given by a polynomial function of the friction parameter obtained by curve 
fitting the results given in Reference 3, Page 175. 
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The flow admittance is also obtained from the results shown in Reference 3. 

ABPqj 

Ad= 

B is a function of (4fL/p) for choked flow and a function of 4fL/p and ^ 02/^01 for 
unchoked flow. The dynamics model uses B in the form of polynomial equations. 


The friction factor, 4f, is a function of the flowrate itself, and hence an iterative 
scheme is required to determine admittance. Starting with an initial guess at 4f, 
the flow Reynolds number is calculated to give a revised estimate. If the flow is 
laminar (Re <2100), the friction factor is given by 


4f = _§! 
Re 


(7) 


If the flow is turbulent, 4f is calculated using a polynomial curve fit of friction factor 
data for smooth pipes. 


Pressure calculations: The pressure change in engine volume I between 
time t and time t + At is given by 


P (I) 

0 t+At 


= P 


»® t 


lit 


C(I)t+ At 


E w (J) 

all J 

connected to I 


( 8 ) 


C (I) is the volume capacity defined as follows 


c ®t + iit 


V (I) 

E T o <V<it 


(9) 


The temperature in volume I when volume I gas is unignited at time t + At is given 
by 


. E q (j) t 


T o® t+ At 




H(I) + At - - „ 

. t all J connected to I 

W(I) t + AtE w(J) t ^3.5 - 3.26 M (I) t+At 


)1 


(10) 


The parameter q (J) is the thermal energy flow through restriction J given by 
q (J) = T q (J) w (J) (3.5 - 3.26 m (J)) 


11 



S-1220 


for cold flow. For ignited flow, q (j) is calculated using a polynomial equation in m(J). 
T 0 (J) and m (J) are the total temperature and percent 0£ associated with restriction J, 
i.e. , those properties existing in the higher pressure volume of the two which J connects. 

It is assumed that only the ignitor and the combustor can sustain ignition. After 
ignition is detected, the temperature and molecular weight in these volumes is determined 
using polynomial equations which were based upon chemical equilibrium calculations. for 
the GO 2 /GH 2 system . 


In order to simulate valve opening time and valve sequencing, the dynamics model 
has provision to express the restriction area as a linear function of time. For time less 
thamthe electrical ON time or greater than the sum of the electrical OFF time and the 
closing time, the restriction area is zero. For time greater than the sum of the electrical 
ON time and the opening time, but less than the electrical OFF time, the restriction area 
is that given in the input form. During opening and closing, the restriction area is assumed 
to be a linear function of time. Figure 3 shows restriction area as a function of time. 

Time increment calculations: The time increment, At, used to advance time in 
the dynamics model, is determined each pass through the model to assure that the solution 
remains stable. The stable time increment is derived below. 


Equation (8)can be rewritten to express the pressure in Volume I at time (t + At) as 
the weighted sum of the pressure in I at time t and the pressures in all volumes connected 
to I 


p o (» 


t+ At 


P o®t 


1 - 


At 

C <VAt 


E 

j 


A d (J) 


At 

C-Wt+At 



where restrictions J connect volume I to volumes K. 


The coefficients of P Q (K) t are always positive, but the coefficient of P ot (I) will be 
negative if 


At > 


£a d <J> 


( 12 ) 


A negative coefficient implies that the higher P Q (I) is at time t, the lower it will be at 
time (t + At ) which is physically absurd. Therefore, At is made small enough to avoid 
negative coefficients in the system by assuring that 


At 


iS 


min 
all I 


C (I) t +At 

T, •'„<•» 

all J connected to I 


(13) 
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Dynamic performance calculations: Thrust, specific impulse, and integrated 
specific impulse are calculated at each increment of time. The thrust is calculated 
using: 


F - P A* C 
o F 


( 14 ) 


Thrust coefficient, Cp, is input or, if combustion has been called, is passed from the 
performance routine. 


Accumulated impulse if given by 


I (t) 


N 


r F.(At). 


(15) 


Integrated specific impulse if given by 

I /f> - 

Sp( ) w* (t) (16) 

Injection analysis. - The injection model provides starting profiles to the 
combustion module. The user must input the injector face diameter, the number of 
injection elements, the propellant temperatures, and the required mixing efficiency. Total 
mass flow, total pressure, and overall mixture ratio are either user supplied or passed 
from the dynamics model. 

The injection model calculates a flow field consisting of concentric annuli of 
uniform flow alternating between fuel rich and oxidizer rich regions. The flow conditions 
are determined to satisfy continuity, conserve energy and momentum, and to achieve a 
desired degree of injector-inducted mixing. 

The injection model analysis begins with the assumption that the sum of the oxygen 
mass flows in the oxygen rich annuli is equal to the sum of the oxygen mass flows in the 
fuel rich annuli. 


w 


O, 


w 


Or 


w 




w 


Hr 


1/2 m w T - 6 


1/2 mw T + 6 
(1 - m) 1/2 w T + 
(1 - m) 1/2 w T - 


_6_ 

m 


m 


( 17 ) 
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The 6 is determined by iteration to yield the input mixing factor defined as in Reference 14 



It is assumed that both oxidizer rich and fuel rich regions have the same total pressure, 
the same ratio of specific heat and have the same static pressure and Mach number. 
Therefore, the ratio of total fuel rich flow area to total oxidizer rich flow area is given by 



The total temperature in the fuel rich region is 


.24 m T + 3. 5 T„ - 3.5 m T 
T ' °2 H 2 H 2 

O — 


( 20 ) 


3.5 - 3.26 m' 


and the molecular weight is 


MW' = 


64.51 


2.016 m' + (32 (1-m')) 


( 21 ) 


Similar expressions apply in the oxidizer rich region. 
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To find velocities in the fuel and oxidizer rich zones, the following mass flux 
expression is solved for Mach number using the Newton-Ralphson method: 


w ' + w 
°2 H 2 - 


( P ° 3 t ') ( 


1 - — 
5 


2 \ - 5 
M \ S' 


) 


( m i. ) 

\ 1 + .2M 2 / 
\ ' 


1/2 


( 22 ) 


where 


P (MW<) 
P = o 


(18540) T ' 


-Jy* V 

f MW' 


up* 

’ Vi + . 


2M 2 


(23) 


Again, similar expressions apply to the oxidizer rich zone. 

The number of flow annuli used in the injection model is determined by the 
following heuristic rule based on the number of injection elements: 


N A 

N e = 1 + ^ 6n 
n = 1 


(24) 


The width of each annulus (and the radius of the inner circular region) is simply the 
injector radius divided by N A . Each annulus contains a region of oxidizer rich flow and 
a region of fuel rich flow. The radii of the discontinuities are determined to satisfy the 
flow area requirement, equation 19. 


After the physical flow field has been determined, the injector model converts to 
the stream function coordinates used in the combustion model. For axisymmetric co- 
ordinates used in the injection model, the von Mises transformation is defined by 




puydy 


(25) 


Once the stream function coordinates of the flow discontinuities have been calculated, 
point flow conditions for the input number of combustion model grid points can be deter- 
mined. The grid points are separated by equal stream function increments. It is 
recommended that, for good resolution, about (10 x N A ) grid points be specified. 
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Ignition. - The ignition model is an integral part of the dynamics calculation 
designed to predict sparking, ignition, and quenching phenomena. If a pilot is present 
in the dynamic system, sparking and ignition is tested for in the pilot and a test is made 
for flame quenching in the combustor. If there is no pilot in the dynamic system, the 
sparking and ignition tests are made on combustor conditions. 

A spark is assumed to occur when voltage is applied to the plug and the potential is 
sufficient for breakdown. Breakdown potential, it is assumed, is a function of spark 
gap and local pressure. Second order effects of gas velocity, electrode design, and gas 
mixture ratio are not considered. Figure 4 shows breakdown potential versus the product 
of pressure and spark gap following Paschen's Law modified by published test results 
(Reference 4). This criterion is utilized in the ignition model as a sparking test using a 
polynomial curve fit to describe the breakdown line. 

If a spark occurs , a test is made to determine if the spark energy is sufficient to 
ignite the ambient gas. The primary factor influencing minimum ignition energy is gas 
mixture ratio. Figure 5 shows minimum ignition energy as a function of mixture ratio. 

This curve, in the form of a polynomial equation, is used as the ignition criterion. Figure 
5 was derived using numerous test results, some of which are shown. Secondary effects 
which are not included in the ignition criterion are gas velocity, electrode design, and 
spark gap. 

If the dynamics model contains a pilot, a test is made to determine if the flame will 
be sustained in the combustor. Figure 6, the result of work reported in Reference 4, is 
used as a quenching criterion. The quenching parameters are mixture ratio, pressure, 
and combustor diameter. The flammability threshold is used in the ignition routine in 
the form of a polynomial equation. The second order effects of flame temperature and 
wall temperature are not considered in the quenching model. 

Mixing and combustion model. - 

Description of techniques: 

1. Describing equations . - The starting point for the mixing and combustion 
model is the boundary layer form of the conservation equations for global mass momentum 
and energy and element and species diffusion. A solution of this system provides the details 
of the flow field including the velocity, temperature, and species fields. 

The global continuity equation can be eliminated from the system of differential 
equations by introducing the von Mises coordinates as the independent variable;. The 
transformation (x,r)-(x, i') is defined according to the relations: 

N N 

Pur = • (26) 

N N 

-pvr = * 
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FIG. 6 
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where N = 


0 - plane two-dimensional flow 


1 - axisymmetric flow 
Introduction of these equations into the differential equations results in: 
Element Conservation 


a/#, ~ [ (Le/pr, (p uW\ 2N ^ ] 


( 27 ) 


Bx 


L. . 


where 


= t VjL oyv a. 


(28) 


and v ^ is the amount of element j in specie i and the W's are the molecular weights. 
‘ Specie Conservation 

(1//) §*[ (Le/Pr) IpUfiAS ) r " r^- j + V° U (29) 

where w i is the volumetric rate of production of specie i. Note also that: 






(30) 


Equations (27) and (29) are used depending upon whether equilibrium or non- 
equilibrium chemistry is considered. If equation (27) is used then i-j equilibrium 
relations are required which are supplied by the indetermancy approached by the 
production term, Wj, as equilibrium is attained. This formulation is used in 
connection with the complete combustion, or diffusion controlled, limit. Equation (29) 
is used when the full H 2 /C> 2 kinetics is considered. 


Momentum 


o u _ 
dx 


/ _N a r . At n, 2n au -j 

(1/4 > iiL (pU|i/4 )x a* j 


(l/pu) 


do 

dx 


(31) 


21 



S-1220 


Energy 


g - a/#) ^ {(pu,/Pr/) r 25J [f| * (Pr-1) ^ 


ou 2 n 


6X 


oa . _ 


+ E h. (Le-1) — i 
i J 

x- 


i / 


(32) 


To supplement these conservation equations, relations among the thermodynamic 
variables are required, viz., 

State : 


P 


P RT IKaVwJ 

i 


(33) 


also 


H = h + 


u 


(34) 


where 


h = L a i ’n i 
i 


(35) 


with* 


h. 

i. 


= h i (T) 


(36) 


In addition, representations for the turbulent transport coefficients fi, Pr and 
Le, are required as well as specification of the chemical system and its associated 
rate constants. 
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With regard to the transport coefficients,, the numerical analysis has been 
structured in such a way as to provide complete generality in evaluation of these 
parameters. That is, they are computed locally and could ultimately be specified 
as functions of the local values of the mean flow variables. At the present time, 
however, the options which have been provided for include only the following: 

For Le and Pr - any non-zero constant value specified by the user 

For p - (a) any non-zero constant value specified by the user 

(b) a modified form of the model due to Hirsch (Reference 5) 
which can be written 


= k 


puD 

1+x/S 


(37) 


where k is a constant which can be input by the user, D is the local thrust chamber 
diameter and S is the injector element spacing which is also a program input. 

For the chemical, system we note first that the Wj are given by 

w ' N p N a . V ”1/ '■ _ 

R • ® N "ip n (_i) ip ^ 

V *i s Mp-VW L J 

P=1 x ~ x 1 


N K 

m = E v'. : N = .E [u" -v' ) 

n U ip P i=l ip 3.P 


i=l 


(38) 


for a chemical system containing N species entering into R elementary reactions 
given by 


N 

E 

i=l 


V M. 

ip i 




N 


v 


i=l 



M. 

1 


P 



( 39 ) 


The present study employs the reaction mechanism given in Table II which 
involves 8 species entering into a total of 17 reactions. 
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TABLE II. H JO CHEMICAL SYSTEM 

M Z 


Reaction No. 


Reaction 

HO_ + H = O + H 6 x 

2 2 2 

H0 2 + H = OH + OH 6 x 

H0 2 + O . = 0 2 + OH lx 

H0 2 + OH = 0 2 + H 2 0 1 x 

H 2 + OH = II + H 2 0 2.19 x 

OH + OH = O + H 2 0 5.75 x 

H + O = H + OH 1.74 x 

0 2 + H = OH + 0 2.29 X 

H 2°2 + ° H = H 2 ° + H °2 1 X 
H 2 0 2 + II = H 2 + H0 2 2.34 x 

H 2°2 + H = 0H + n 2 ° 3.18 x 

O+H+M=0H+M lx 

0 + 0+ M = 0 '+M- 9.38 x 

H + H + M = H 2 +M 5 x 

H + OH + M = H 2 0 + M lx 

0 2 + H + M = H0 2 + M 1.59 x 

OH + OH + M = H 2 0 2 + M 8.4 x 


Notes : 


k £ =AT B e - E/RT 


“b = V k o 

k is calculated internally by the computer at each temperature utilizing 
C free energy of formation of each reactant 
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To complete this formulation, initial and boundary conditions must be specified. 

• Initial conditions - The initial conditions must represent the details of the 
flow emerging from the near region and, therefore, must allow for the specification of 
velocity, temperature, and composition profiles. Modeling of the initial conditions 
is carried out as follows: referring to Figure 7, the actual pattern emerging from the 
near region is divided into, an arbitrary number of annuli. To remove any three- 
dimensionality, the pertinent variables are circumferentially mass-averaged in each 
annulus. The mass averaging ensures conservation of mass, energy, and momentum. 
Depending upon the particular pattern, smooth rather than stepped profiles may be 
appropriate. 


These conditions are given by: 


@ x = 0 (starting station) and 0 r £ r (0) 


0 s: r 


u = u 1 (r) 

^ T = T x (r) 

= a i,l (r) 


r l< r * r 2 




u = u 2 (r) 
T = T 2 (r) 


= a i,2 <r) 


r k-l <r£r k 


u = Uj^r) 

■T = T k (r) 

a i = a i.k (r) 


• Boundary conditions - th rust chamber geometry - To provide the versatility of 
either predicting behavior of existing hardware of designing new hardware, the analysis 
was developed to permit specification of the chamber contour or the axial pressure 
distribution. Specifying one renders the other a dependent variable. Thus, for x^ 0 


if p (x) is prescribed 

r (x) is computed 
w 

if r w ( x ) * s prescribed 
p(x) is computed. 
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This option provides the capability to evaluate the effect of acceleration or deceleration 
of the chamber flow upon the mixing rate, ignition and flame propagation rate, and 
local wall heat transfer for a given injection configuration. 


• Boundary Conditions - Axis 


for x s o and r = 0 


3a . 

i 

3^ 



3>& o 'I' a 'J 1 



• Wall boundary conditions - The interaction of the chamber flow with the 
wall involves drag, heat transfer, and, in general, mass transfer - the latter being a 
consideration for transpiration and/or film cooled chambers. The various combinations 
of boundary conditions which have been implemented in this computer program are 
indicated in Table HI. Note that for the transpiration cooling model an explicit boundary 
condition for a "tracer" specie is included. This is needed for proper implementation of 
the boundary condition in this case as discussed with the chamber cooling models. The 
manner in which "bulk" values of the several parameters indicated by subscript b are 
evaluated, is also described in the chamber cooling section. 


2. Computational procedures. - The solution of the above system of equations 

is obtained employing an explicit finite difference techniaue. Figure 8 shows a generic 
point, (n+1, m) in the x-\p grid network. The finite difference formulation for the calcula- 
tion of the flow at the point (n+1, m) is obtained by using the following explicit/difference 
relations where P is any one of the pertinent variables. 

dp — J: — (p _ p ) 

3x Ax n+l,m n,m (40) 


££ = — ^—(P _ - P , ) 

d-\jr 2 a\E’’ n,m+l n,m-i 


(41) 


£-b ££ 
b'i oW 


A^' [ b n,n^ (P n,^l b n,ia-^ n,m n,m-l 3 



(42) 
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where 


, . /t N, 2N 

b = (pujj,) /w ) r 

b , , = ~ (b +b ) 

n,m + i 2 2 n,m rwinjl 


and 


= mA'J' 


(43) 

(44) 


(45) 


The conservation equations in difference form are: 

Elements : 
m - 0; 


^j^n+1,0 {a j*n,0 


. . . + 


L (Lep/Pr) (pU) Jn,o[ (a j ) n,l' (a j ) n,0 ] 


(46A) 


m^O 
(a^) 


= (a.) + 


A.X 


j'n+l.tn '“j'n.m- ,N 2»s { (t * b/W n,«HJs S’ «>,«+}. 


\ __ 

... - j (Le b/Pr) ^ + (Le b/Pr) , I (a.) 

L. n, n,m~'2 J yn. 


m 


. . . + (Le b/Pr) , 
. . n / 



L 


n, mr-1, J. 


(46B) 


Note: The species conservation equations have the identical form with the production 
term added to the right hand side. 
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Mnmp.ntnm 

m = 0 


. 2/vx(N+l) r 

u =s n + — • 

n+1,0 n,0 (A'M 


j M (P u ) 1_N " 1 L 0 f u n -u 1 ... 

i / v? H/X n,Oj 


_ - Ax .do 

(pu) Wrx+i 


(47) 


m 7* 0 


u 


= u T 


Ax 


{ b 

L n,i 


n+l,m “n,m N. 2 +n \^n,m+^ U n,m+1 
m (a v K 


...-u (b , + b ,)+b i u , 1 
n,m n,m +% n,m =i 5 11 , 111 -% n,m-l j 


... - ■<**> (^> , 

pu n,m dx n+1 


(48) 


Energy 


m = 0 


1-N_ 


H n + 1.0= H n.0 + >».o] ( <fc> n. 0 <V I'Vo 1 


I (1 -57>n.O ‘»i.l - U n.0> ••• 

- <*i>n,0 ] (h i 17 ] O ,0 }' 


(49) 
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m / 0 


H , . =H + 
n+l,m n,m 


N. a ^{ tb/Pr> .,^ H n,» + l 


m (A 'I') 
-(b/Pr) 


, + (b/Pr) : H 
n,m+H n,m-% n,m 


+ (b/Pr) 


n,m-h "n/rn-l ‘ 2 


«-• — (b-b/Pr) u : 


n , m + > 2 n , m+1 


- 7 I" (b-b/Pr) ■ + (b-b/Pr) , ] uj 
2 L n, m-T'i r^m-^ _j n. 


m 


+ j (b-b/Pr) 


E (a.) 


+ E (a i ) 
i 


1 

u a 


n,ra- 

r (bh. 

Le-1 

L i 

Pr 

-l (bh i 

Le-1 

Pr 


+ E (a,) 




i_ i v Pr Jn,m+V 


)_ + (bh, 


i Pr n,m-^ J 


-h ; * * * 




(50) 


The boundary conditions are expressed in finite difference form by using 
Equation (41). The use of such a central differencing scheme at the wall is implemented 
by carrying along an additional streamline above the wall. 

• Step size control - The step size in the explicit finite difference scheme is 
controlled by a stability criterion and from studies of linear parabolic partial differ- 
ential equations there results the following condition: 


AJ? J" Pr/Le u l 1 rr N a\t, 2+N ~ 

6(1+N) L - 5n '°‘ i 3 fre b/?r) + (Le b/Pr) , 
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Description of cooling techniques 
Slot cooling 

Describing equations: The mathematical model selected for implementation within 
the present program utilizes the Hatch-Pappel (Reference 6) correlation together with the 
Bartz (Reference 7) method for evaluation of the requisite heat transfer coefficient and 
the Sellers (Reference 8) procedure for including the effect of multiple slots. The manner 
in which the latter is accounted for is indicated schematically in Figure 9 where the various 
parameters pertinent to the slot cooling problem are defined. According to the analysis 
cited above, the wall temperature distribution between the Ith and the Ith+1 slot can be 
determined from* 


X I I ITT (51) 

m n = (Pr^T) © [.04 - (x-xp/G s ] 


where the effectiveness parameter is defined by 


and 






(T 1 - 1 - TVer 1 - 1 .- T 1 ) 

w w w c 

(OUC )W 

p c 

1/8 " r 

[1+0.4 arctan (U^-l). ] 

13/8 I [1-1/U 1 ] 




; U a 1.0 
; ' U 1 £ 1.0 


(52) 


(53) 


(54) 


*More precisely, these equations apply in the region bounded by (x-x 1) /G* s* s . 04 
and x ^ x^ + 1 . For axial locations upstream of this lower limit but beyond x^ the 
temperature is constant and equal to T* as has been indicated in Figure 9. 
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Following Bartz, the heat transfer coefficient hg , appearing in Equation (53), 
is calculated from 


where 



1 /A? 

a 


0.9 


- .025 
_ (D* ) 


.2 

c p 

tv 




.8 






(55) 


(56) 


In these relations, the subscript b indicates the bulk properties of the main thrust 
chamber. The manner in which these are evaluated is indicated in the next section. 


For a single slot (or for the first slot) it is necessary to define more explicitly 
the significance of . Consistent with the manner in which the Hatch -Pappel model 
was developed for single slots, we take equal to the adiabatic wall temperature 
based on bulk properties (T aw In turn, this is related, in the manner which will be 
indicated in the next section, to (H aw )b* 

1/3 2 . 

Where ^aw’b = V + (Pr> b V 2 <57) 


Computational procedures: The slot height s*, axial location x* , unit mass 
flux rate of coolant m* , coolant reservoir temperature T J , slot exit pressure p£, 
coolant Prandtl number Pr c and specific heat Cp c is considered to have been prescribed 
for all I so that all parameters with subscript c appearing in Equations (51) through (54) 
have been evaluated. At a generic point x = x n such that x^ £ x n £ x^ +1 it is further 
assumed that all flow parameters throughout the thrust chamber have been determined. 
The problem is to determine the wall temperature at the next integration step x = x n +i 
which will be denoted by (T w ) . The computational procedure utilized to accomplish 
this (assuming x n+1 < xj, + 1 ) is as follows : 

(a) compute bulk properties of conservation variables at x = x^ from 

r 

27r J 0 WpuHrdr 

^ ih 
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U L 




m 


i. 

2 irf W pu a ,rdr 


where r w and m denote the values of thrust chamber radius and total mass flux through 
the motor at x = x , respectively. 

(b) compute all other bulk thermodynamic properties such as p^ p ,, 

(Taw) etc * f rom appropriate auxiliary thermodynamic relations, 

e.g. : 


(T ), = 

aw b 


:cn 


[ (H ). 
aw o 




where (H follows from Equation (57) and the functional notation implies the use of 
the internal enthalpy -temperature fits incorporated by the program.* 


(c) determine the wall composition at x = x n+ j , (a iw ) n+1 , by application 

of the film cooling boundary condition for species diffusion (c.f. Table III). 


Since 



W T, 
w b 


W. 


T 1 


and 

P\v ~~ ^ c n (T w ) 

the only unknows appearing in Equation (51) are (T w ) n+ ^ and (T ) 


*It is important to note here that p 5 represents a laminar (molecular) viscosity which 
depehds 6 n both the composition and temperature. Although more complicated representa 
tions for the viscosity of a mixture could be incorporated to evaluate this parameter, it 
has been deemed sufficiently accurate for the present purpose to ignore the dependence on 
composition and evaluate p from the Sutherland formula for air; viz: 

- , , , .. ^-8 ib.sec. 

** b 

where is in degrees Rankine. 


= 2.27 


T b 3/2 

T.+1S8 

b 


x 1 C' 
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But 


1-1 

•tn 7} 


1-1 1/8 
(Pr R ) 
c c 


o 1 ,;L r .04-(x-x_ I " 1 )/G I “ 1 s I 1 ] 


1-1 

V 



1-2 


1-1 

•T 

w 


)/(T 


1-2 


w 


■T I " 1 ) 
c 


etc. , where the subscript n+1 has been omitted for clarity. Thus, the next step in the 
procedure is to 


(d) systematically determine in order 


(T ) , 

v w' n+1 


(T ) 2 

' w' n+1 


<T ) I_ ] 
w ; n+1 




by iterative solution of the transcendental Equation (51). 


Solution of the last member in this sequence gives the desired value of 


(T ) 
w 


I 

n+1 * 


Transpiration cooling 

Describing equations: The mathematical model selected for implementation within 
the present program uses the Bartle-Leadon (Reference 9) correlation together with a 
reference state approach to evaluate the requisite Stanton number. The reference state 
is defined by utilizing the Eckert (Reference 10) definition of reference enthalpy and 
Knuth's (Reference 11) representation for the reference composition. In accordance with 
these formulations, the temperature of the transpired surface is related to the coolant 
flow rate by the relation: 


where 


(T 

(T 


aw o 


aw o 


T 

w 

T 

c 


i - u + §) 


G 



(me ) 
p c 

(ouc ) 

• P e 


(58) 


(59) 
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and 


(T ) 

aw o 


(H ) 
aw o 


fcn L , (H aw } o ' a ie ] 
he + (Pr) 1/3 -f~ 


Here the subscript o indicates properties evaluated under "no blowing" conditions, while 
subscript e indicates edge conditions. Evaluation of the latter is discussed later. 


Use of a Reynolds analogy and reference state in conjunction with an incompressible 
skin friction law leads to the following representation for the Stanton number: 



-2/3 

.0592 (Pr) • (£_) 

e pa 


p*u x 



-i/5 


(60) 


where the asterisk implies that the thermodynamic variables are to be evaluated at the 
reference state condition. The latter is defined by combining Eckert's form for the 
reference enthalpy 


h* 


“ (H +h ) + 0.22 (H -h ) 
2 w a e 


(61) 


with Knuth's representation for the reference composition: 


a* 

x 


-v X 

■ “i 

e c 


i / c 


where 



In (w /w ) 
e w 


l (i--a ) w / u- a ) vV ] 

c e c w 

e w 


(62) 

(63) 


W = (1-a ) - (i: 

e c . / 

e x^c 


a./WV) 
x x 


-1 

e 


W 

w 


(E a./w. ) " 1 

X X w 
X 


We note here that the parameter W g represents the molecular weight of the subsystem 
consisting of all molecular species except the injected species (coolant) at the edge of 
the boundary layer. 
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Computational procedure: The coolant reservoir temperature T c (x) is considered 
to be given and either the coolant mass flux rate (x) (Option A) or the wall temperature 
T (x) (Option B) are prescribed. Assuming that all flow properties throughout the thrust 
chamber are known at the generic point x =x n » the problem is to determine either T at 
x = x n + i or m at x = x n +1 . From the point of view of integrating the describing 
equations, these two options differ in a very fundamental way. Thus, for Option A (m c 
specified; T to be determined), Equation (58) is directly coupled to the system of 
finite difference equations which are to be integrated. In contrast, for Option B 
(T w specified; m c to be determined) the integration could proceed independently of 
whether or not a solution of Equation (58) is obtained. The procedure used in the case 
of Option A will be described first. 

Option A 

(a) Compute "edge" properties at x =x n * Provision has been made for 
evaluation of these edge properties in two different ways. One of 
these is identical to the procedure used for the film cooling model 
wherein the edge conditions are equated to the bulk properties of the 
conservation variables. The alternate procedure involves determination 
of the edge of the boundary layer by establishing the extent to which the 
tracer gas has penetrated laterally.* At this point in the flow field, 

the values of H, u and cq can be determined from which all parameters 
appearing in Equations (58) through (63) with subscript e, can be 
- evaluated. 

(b) Determine the wall composition at x = x n+ i, («i w ) n+1 , by application of 
the transpired wall boundary condition for species diffusion (c.f. 

Table III). Since m c at x = x n is known, this can be readily accomplished. 

(c) Determine the reference composition * at x = x n+1 , from Equations (62) 
and (63). Note here that in this determination we are "lagging" on the 
edge condition. 

(d) Determine T w at x = x n+1 , by iterative solution of Equation (58). Note that 
this is required since the parameter G depends implicitly upon T w , via 
Equations (59), (60) and (61) and the fact that H w is a function of oq and 
T w at x = x n+1 • With T w and, therefore, H w determined at x = x n+1 , all 
data needed to continue the integration of the conservation equations are 
available. 

Option B 

For Option B, steps (a) through (c) are identical. However, since T w is known 
at x = x n +i , at that point can be evaluated directly which allows h* and 
therefore p* and p* to be evaluated, using Equation (61). The only unknown 
appearing in Equation (58) therefore is the mass flux rate at x = x n+1 which 
can be solved for in a straightforward manner, if desired. 


*The criterion used to determine this penetration is to locate the point wh6re <Xj, s . 01 ttT w * 
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Performance model. - The output of the combustion model provides the charac- 
teristics of the gases at the throat of the rocket engine in terms of mixture ratio and 
enthalpy along each streamline. The data thus obtained can be used to determine the 
specific impulse and thrust of the engine as a function of the nozzle characteristics. The 
output of the combustion model can be used as input to the various CPIA standard programs 
(Reference 1), but because of the complexity and relatively long run time of these pro- 
grams, a more simplified method was used to determine the losses in the nozzle. The 
losses are due to: 

(1) Kinetic recombination 

(2) Boundary layer or viscous effects 

(3) Divergence 

(4) Film and dump cooling 

The performance of the rocket engine can be estimated using approximate methods 
as described in Reference 1 for the boundary layer or viscous losses and the divergence 
losses, while the kinetic recombination losses are determined for each streamline using 
the methodology outlined in Reference 12. The enthalpy and mixture ratio at each stream- 
line is used as the input to determine the specific impulse after kinetic loss. The char- 
acteristics of the streamline are used in determination of the viscous loss, while the 
loss due to divergence is based solely on the area ratio and an infinite pressure ratio 
(y = 1.20). The losses due to kinetics and mixture ratio are thus mass averaged and 
the viscous and divergence losses subtracted directly. In addition, any hydrogen or oxygen 
dumped in the expansion bell is assumed not to mix and the I of the cold gas is assumed 

bp 

to apply. In summary, the equation used to obtain the performance is 



Kinetic recombination: The method used to determine the amount of recombination 
which takes place in the nozzle for each streamline is based on an approximate method 
developed by United Aircraft Research Laboratores and is described in Reference 12. The 
losses are evaluated using the Modified Bray criterion for predicting the point in the re- 
combining nozzle flow where the reactions have departed significantly from equilibrium. 
This criteria has been successfully used to analyze flow in which only one reaction is 
kinetically important. The Bray method was extended to multicomponent-multireaction 
performance calculations and the method was incorporated in the performance model. The 
basic method use d^to^e ter mine the kinetic lossddr _ each=streamline is - as described below. 

(a) The specific impulse for each streamline based on total pressure, enthalpy, 
and mixture ratio in both equilibrium and frozen flow is determined. 

(b) Based on the nozzle contour and design, the amount of recombination 
at mixture ratios of 4,6, and 8 are determined at the known pressure. 
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Assuming that the frozen and equilibrium specific impulse below a mixture 
ratio of 2 are identical and that the percent of recombination above O/F = 8 
is the same as O/F = 8, the data for the recombination in (b) was linearly 
interpolated for each streamline. 


The data presented in Reference 12 was based on liquid O 2 an d H 2 . In order to 
make the data adaptable to propellant at any temperature, the data was modified to the 
form of percent recombination of 

Is ">K - H 

’’R * 

■sPs ' W F 

Thus, if the equilibrium and frozen specific impulse at a known enthalpy and mixture ratio 
are known, the kinetic data can be readily determined. 


Figures 10, 11, and 12 show typical data used in the determination of the kinetic Igp 
The nozzle design area ratio is obtained from the Figure 10 based on the dimensional 
characteristics. The nozzle gradient (See Figure 11) is then determined for various area 
ratios (in the case of this analysis from 1.01 to 5.0). The gradient value is then divided 
by the nozzle throat diameter and multiplied by the line labeled "Equilibrium" obtained 
from Figure 12. This data of H atom gradient vs. area ratio is then cross-plotted against 
curve labeled "Composite Kinetic (Transition Factor = 1)" of Figure 12. The intersection 


is the freeze area ratio. In the reference, the kinetic I sp is found directly from a plot of 
I g p vs. freezing area ratio (freeze) (see Figure 13), but as described previously, the 
data was modified to yield percent recombination. 


Boundary layer or viscous losses: The losses due to the interaction between the 
gas and the wall, commonly called boundary layer losses, are determined using the 
approximate methods described in Reference 1. This method is based on the results of 
a more rigorous calculation method, commonly called TBL (turbulent boundary layer). 

Data for the various characteristics of the boundary layer, momentum thickness, for 
example, have been presented in terms of the isentropic exponent, the temperature ratio 
(gas-to-wall), and Mach number. In the present analysis, the isentropic exponent, 7 , 
was assumed to be 1.2 and the temperature ratio 0.2, since a majority of the wall tempera- 
tures are in the 500-2000°F range while the gas temperature is >6,000°F. Figures 14, 

15, 16, and 17 show the data used. The effect due to variations in contraction ratio and 
radius of curvature was assumed constant since the effect is second order. The equation 
thus used to determine the loss in therms of thrust is 
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FIG. 10 
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(REF 12) 
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FIG. 12 
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NORMALIZED GRAPHICAL SOLUTION FOR FREEZING AREA RATIO 
USING MODIFIED BRAY ANALYSIS (REF 12) 
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FIG. 17 
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where 



( 66 ) 


Dividing by the total mass flow results in I g p loss. The gas properties were taken to be 
those of the streamline adjacent to the wall. 


Divergence: The loss due to divergence was obtained from Reference 1 and is 
shown in Figure 18. The pressure ratio was assummed to be infinite. 

Mixture ratio maldistribution effects: The gas properties at each grid point from 
the combustion program must be integrated to find an average kinetic specific impulse 
as described by the first portion of equation 64. The for each streamline is a point 
function which is applied to a flow streamtube. The method used is simply one of linear 
interpolation between streamlines. For example, the following data 

Streamline Radius Area (i) of Streamtube 
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FIG. 18 


DIVERGENCE EFFICIENCY FOR CONTOURED NOZZLES 


P/P =0 
e o 

■y = 1.2 

from Ref. 1 
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where r n = wall radius, will, when combined with the individual ( w/A)^ from the 
combustion program provide the \v i required for equation (64). 

Heat transfer model . - The heat transfer model predicts combustion, injector, 
and coolant temperatures for a generalized G0 2 /GH 2 rocket engine. The cooling methods 
which can be used are shown in Figure 19. The user specifies combustor geometry and 
heat transfer properties, engine operating conditions, and, where applicable, coolant 
injection conditions. The heat transfer subroutine creates from the input a thermal model 
of the engine and solves for temperature distribution using a finite difference technique. 

The heat transfer model is completely unconnected to other characterization program models. 
It exchanges no information with other models. Unlike the other characterization sub- 
routines, heat transfer input is accepted in the subroutine itself, not passed from the driver 
program. 

The thermal model created by the heat transfer routine consists of mass nodes 
connected by thermal resistances analogous to an electrical network containing capacitors 
and resistors. A simple example is given in Figure 20 showing a bar heated by convection, 
cooled by radiation, and exhibiting internal conduction. It could represent a portion of 
film cooled combustor wall. Figure 20 also shows how the thermal admittances (admittance 
- (1 ~ resistance)) are calculated for conductive, convective, and radiative heat transfer 
and how the thermal capacities are defined for mass nodes and for boundary nodes. Surface 
areas, cross-sectional areas, lengths, volumes, view factors, and convective coefficients 
are automatically computed by the heat transfer sub-program in the process of creating 
the thermal model. 

Heat transfer solution technique: The basic equation describing the continuity of 
heat flow at node i. is 



connected 
to i 



T i 

t + At 

At 



which can be rearranged to give 


T i 


t + At 


explicitly. 


(67) 


T i 

t + At 



connected 
to i 


( 68 ) 
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FIG. 19 


HEAT TRANSFER MODEL 
COOLING METHODS 
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FIG. 20 


HEAT TRANSFER MODELING 
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In other words, the temperature of node i at time (t + A t) depends upon the weighted 
average of the temperatures of nearby nodes and the temperature of i at time t. 


Except for the coefficient of , the weights are inherently positive. In order to 
assure solution stability, all Tj coefficients are kept non-negative by selecting time 
increments, A t, such that 


C- 

At ^ Za~ f ° r aU 1 (69) 

all j ^ 
connected 
to i 


For steady state calculations, the temperature at node i is simply: 


T. =E 
1 



£ A 


D 


i> J 


T 

J 


all j connected to i 


(70) 


Calculation of passive cooling heat transfer admittances: Hot side convective 
coefficients are calculated using the method of Bartz (Reference 7) shown in Equation (55). 


Viscosity is calculated using Sutherland's formula 


^ (46.6x10 ) y MW (T 0 ) lbm/in sec 


.6 


Prantl number is calculated using the approximate result 


Pr - 


4 y 


(71) 


97-5 


Specific heat, Cp, is utilized in the form of a polynomial curve fit of rocket 
performance data evaluated at a mean film temperature.. 


Radiation view factors are assumed to be 1.0 on the external surfaces of the 
engine. The inside of the combustor is assumed to have a negligible view of the engine 
surroundings. For the purpose of calculating radiation view factors, the inside surface 
of the expansion bell is assumed to have approximately the shape of a 15® half angle cone- 
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The local view factor is given by a polynomial curve fit of previously calculated view 
factor data for 15° cones. 


y = 5 < l/ d* . x a> 


Active cooling models 


Regenerative cooling: The coolant side convective heat transfer coefficient is 
calculated using the equation developed in Reference 13. It is assumed that only hydrogen 
is used as a coolant. The coefficient is given by 



(72) 


whe re : 


Re b 


. D 
w H 

A P 


Hydrogen conductivity and viscosity are determined using polynomial curve fits 
of Reference 13 results with bulk coolant temperature as the independent variable. 

Coolant temperature rise is calculated using heat transfer resistances connecting 
the fluid nodes and given by Res = 1/wCn- This "heat flow" model is modified to prevent 
upstream node temperatures from being affected by downstream nodes, much like a 
cathode follower in electronics. The coolant temperature rise, then, is calculated using 
the same modeling framework shown in Figure 20. 

Film cooling: The film cooling model is based upon experimental and analytical 
results given in Reference 6 for a variety of injection conditions and described in the section 
on the combustion model. 


Coolant velocity is determined using input parameters and mass continuity. Coi’e 
gas velocity is calculated from core Mach number which is determined by solving 



1_ 

M 



i 7 + 1 

2(7 - 1 ) 


(73) 


for Mach number by Newton's method. The ratio of specific heats, 7 , is available as a 
polynomial equation based on equilibrium chemistry results. 


The film cooling model is used for multiple injection cases by using the film 
temperature produced by upstream injection as driving potential for films injected down- 
stream. Calculations for the film injected farthest upstream are carried to the nozzle 
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exit using core temperature as the driving potential. Calculations for the film injected 
next downstream are based upon the film temperatures resulting from the first film. 
Temperature distribution for a third film is influenced by the two upstream films, and 
so forth. 


Liner cooling: The liner cooling model is designed to predict coolant temperature 
rise, and combustor temperature distribution. The liner is. assumed to be a thin metal 
structure which delivers hydrogen coolant downstream. It is cooled by convection and 
radiation on the backside, and, in addition, may be film cooled on the combustion gas 
side. 

Backside convective coefficients are calculated using equation (72). Since liners 
are sometimes run quite hot, radiative heat transfer to the combustor is important. The 
effective emissivity between liner and combustor is given by: 


C = 


1 

^ liner 


1 


+ 


1 

com bus tor 


1 


(74) 


which accounts for reflection back to the liner. 


Like the regenerative coolant temperature rise model, the backside liner flow 
model is modified to allow "heat flow" only in the downstream direction. The liner efflu- 
ent forms a cooling film for the combustor downstream of the liner. This film is treated 
exactly as described in the film cooling section. 


Combined cooling models: The basic cooling models may be combined in a number 
of ways to describe complicated cooling schemes. The radiation and chamber conduction 
heat transfer modes are always in effect unless deactivated by setting emissivity or con- 
ductivity equal to zero. Film cooling may be used in conjunction with regenerative cooling 
or liner cooling. There is no provision for combination liner and regenerative cooling 
nor for transient calculations with active cooling. 


Injector heat transfer model: The program contains an option for including a 
thermal model of the injector. Due to the vast number of possible injector configurations, 
most heat transfer information is required as input including heat transfer coefficients, 
surface areas, face heat flux, injector weight, injector-combustor seal resistance, etc. 
The model considers heat input to the injector from face convection, conduction from the 
combustor, and radiation from the combustor. Cooling is by radiation and by convection to 
the flowing propellants. The injector option may be used in both steady state and transient 
operation and in conjunction with any of the active cooling models. 
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SAMPLE CASE RESULTS 


The characterization program was tested using an engine for which test data was 
currently available. This engine was fabricated and tested by the Aerojet Corporation and 
the results were presented in Reference 2. The results of the test case indicated that 

a. The performance and dynamic response of the engine could be 
simulated by the computer program 

b. The heat transfer analysis was adequate for predicting wall and coolant 
temperatures.. 


Dynamics and Ignition 

The Aerojet engine used in the dynamic analysis is defined in Figures 21-23. 

The initial conditions and the calculated volumes and valve response are shown in 
Figure 2. The sample case output is shown in Appendix B. The plotted results are shown 
in Figure 24. The response of the actual engine and the output from the computer simulation 
are essentially identical. 


Injection, Combustion, and Peformance 

The results of the performance calculations which use the output of the combustion 
and injection model indicate that the methods used to analyze performance can predict the 
specific impulse within the accuracy required. The results of the study are shown in 
Figure 25 where data from Reference 2 is plotted against the output of the sample case. 


Heat Transfer 

The heat transfer program was used in a comparison of the engine cooling scheme 
shown in Figure 26. The engine which is both regenerative and slot film cooled has not 
been tested and as a result, the comparison could only be with the analytical results re- 
ported in Reference 2, shown in Figure 27. As shown in Figure 28, the results of the 
heat transfer computer program are essentially compatible with the results of Reference 2. 
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IMPINGING ELEMENT INJECTOR MANIFOLDING 
AND FLOW SCHEMATIC 
. (REF 2) 


FIG. 21 
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FIG. 22 


HEAT SINK CHAMBER 
L* - 15 INCHES 
(REF 2) 
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FIG. 23 
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FIG. 26 


AEROJET RE GENERA TIVELY COOLED 

THRUST CHAMBER 
(REF 2) 

(Used for Heat Transfer Analysis) 















CONCLUSIONS AND RECOMMENDATIONS 


A FORTRAN V computer program was written which predicts the performance of 
G0 2 /GH 2 rocket engines during both pulsing and continuous operation. The results of 
the comparison between the test case used for the computer simulation and the actual test 
data indicate the validity of the methods used. The adequacy of the comparison is de- 
pendent upon the input used. Critical parameters such as injector mixing efficiency and 
turbulent viscosity (used in the combustion program) must be chosen with care. 

It is recommended (1) The program be tested with at least 3-4 engine concepts for 
which data is now available, and (2) The injection model be modified to include the data and 
models now being generated under NASA Contract NAS 3-14379 (Investigation of GH 2 /G0 2 
Combustion). 
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SYMBOLS 

speed of sound (total conditions) 
flow area (geometric) 
thrust chamber throat area 

mass flow admittance or heat transfer admittance 

function of (4fL/D) 

skin friction coefficient 

thrust coefficient 

specific heat 

wc of Ith node 
P 

pA */m 

orifice coefficient 

thrust chamber throat diameter 

diameter 

Fanning friction coefficient 
thrust 

mass transfer parameter for Ith slot 

mass transfer parameter for transpiration 

static enthalpy of mixture 

static enthalpy of Ith species 

heat transfer coefficient for Ith slot 

reference state enthalpy 

total enthalpy 

empirical constant 
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I 

Isp 

k 

Vp 

k 

c,p 

k 

f,P 

L 

Le 

M 

m 

c 

m 

MW 


m 


E 


N 

n a 

P 

P o 

Pr 



r 

r 

w 


SYMBOLS (Continued) 

impulse 

specific impulse 
thermal conductivity 

backward, equilibrium, and forward rate constants 
for the pth reaction 

length 

Lewis number 

Mach number 

coolant mass flow rate 

thrust chamber mass flow rate 

molecular weight 

percent oxygen 

number of injection elements 

number of annuli 

pressure 

stagnation pressure 
Prandtl number 

heat transfer rate from thrust chamber to wall 

radial coordinate 

local thrust chamber radius 


70 



S-1220 



c 


R 

R 

R* 


Re 

Res 

I 

s 

St 

t 

T 

T 

o 

u 


V 

U 


I 


w 


w. 


\v 


w. 

1 


W 


X 


SYMBOLS (Continued) 

( ous/jjl) * Reynolds number for Ith slot based on coolant con- 
ditions and slot height 

gas constant 

universal gas constant 

thrust chamber throat radius 

Reynolds number 

resistance to heat flow 

height of Ith slot 

Stanton number 

time 

temperature 
total temperature 
axial velocity 
volume 

(V 

mass flow rate 

volumetric rate of production of species i 

total mass flow 

molecular weight of species i 

weight 

axial coordinate 
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x. 


a 


a 


a * 

1 


A F 

I 

T? 


BL 


p* 




M-w 

P 

Oo 

p* 

a 1 

y • 


SYMBOLS (Continued) 


pxial coordinate of Ith slot 

mass fraction of ith species 

mass fraction of jth element 

reference state composition 

boundary layer loss 

effectiveness parameter for Ith slot 
divergence effectiveness 


turbulent viscosity 

molecular viscosity based on reference state properties 
configuration factor 

molecular viscosity based on bulk properties 


molecular viscosity based on wall conditions 
density 

stagnation density 

density based on reference state properties 
thermodynamic parameter 
stream function 
specific heat ratio 
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SYMBOLS (Continued) 


Subscripts 


BL 

boundary layer 

e 

edge conditions 

w 

wall conditions 

b 

bulk conditions 

H 

hydraulic 

aw 

adiabatic wall 

c 

coolant 

T 

trace 

o 

zero blowing condition 

D 

duct 

V 

volume identification (typ) 

At 

time increment 

K 

kinetic 

S 

shifting equilibrium , 

F 

frozen equilibrium 


Superscripts 

' fuel rich region 

" oxidizer rich region 

— average 


refers to throat 
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S-1220 


O /H characterization program input parameters 

£ Z 


NAME 


MODEL 


FORMA T 


DESCRIPTION 


NDYNA 


Driver 


NCOMBU 


NINJE 

NHEAT 


I 3 


1 if dynamics calculation, 0 otherwise 
1 if combustion calculation, 0 otherwise 
1 if injection model used, 0 otherwise 
1 if heat transfer calculation, 0 otherwise 


NNODE 

NCONN 

NITER 

LPRTF 

NCHEK 

PCHEK 

PSST 

NCOMCL 

DPCOMC 

CD 


Dynamics 

tf 
• H 
H 

ft 

H 

! f 

M 

1! 

tf 


I 4 


t! 

E10.2 


ft 


No. of nodes (volumes) 

No. of connectors (orifices or ducts) 
Maximum no. of iterations 


Print frequency (print results every LPRTF itera- 
tions) 


Time is advanced to the first valve closing time when: 


dfpRES(NCHEK)) 
d (TIME ) 


< PCHEK (psia/ms) 


once PRES(NCHEK) a PSST (psia) 


I 4 

F10.2 


1 1 


Call combustion every NCOMCL dynamics iterations 

Call combustion every time combustion pressure 
increases by DPCOMC (psia) * 

can use either or 

Nozzle coefficient both 


NOINJ(I) 

NVOL(I) 

ICOMB(I) 

VOL(I) 

PRES(I) 


24 I 3 

I 3 


F10.0 

F10.0 


Volume numbers feeding injector face (used by 
injection model) 

Volume number (consecutive integers , 1 = 1,...., 
NNODE) 

ICOMB(I) = 0 Neither pilot nor combustor 
" ' 1 Unlit combustor 

2 Lit combustor 

3 Unlit pilot 

4 Lit pilot 
3 

Node volume (in ) 

Initial pressure in Node I (psia) 
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O /H CHARACTERIZATION PROGRAM INPUT PARAMETERS 

— 2 *- - 2 ■ 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

TEMP(I) 

Dynamics 

F10.0 

Initial temperature in Node I (°R) 

RMK(I) 

ft 

If 

Initial Massfraction O in Node I 

IADMIT(J^) 

tl 

I 3 

Volume no. connected to admittance J (I = 1,2) 

IRTYPE(J) 

ft 

ft 

IRTYPE(J) = 0 orifice 

= 1 duct 

= 2 throat 

CAREA(J) 

1! 

F10.0 

2 

Full open cross sectional area (in ) 

CCOEFF(J) 

ft 

f t 

Discharge coefficient for I = 

0, 2 

DLEN(J) 

f t 

ft 

Length (inches) for 1=1 


TIMON(J) 

ft 

' 

Time valve starts to open(sec) 

If all timing parame- 

TIMOFF(J) 

ft. 


Time valve starts to close(sec) 

ters are 0.0, area is 
• assumed constant at 

TIMOPN(J) 

ft 


Valve opening response time 

CAREA 

TIMCLS(J) 

tf . 


Valve closing response time 


SPTIME 

ft 

It 

Time of initial spark (sec) 


SPKTL 

T» 

ft 

Time of last spark (sec) 


SPKF 

If 

It 

Time between sparks (sec) 


SPGAP 

ft 

ft 

Spark plug gap (inches) 


SPARKP 

If 

ft 

Spark plug potential (volts) 


SPARKE 

II 

ft 

Spark energy (millijoules) 


DC 

ft 

It 

Chamber diameter (inches) 


- 

Combustion 

12A6 

Title card - any statement. Will be printed on each 
radial profile output. 

MPSI 

IT 

I 5 

Number of grid points at x = 

x initial . 

IPRESS 

ft 

M 

Number of grid points at x = 0, used in halving 
the grid, and program restarts at x >0. 

ITURB 

ft 

ft 

1 - Use Hersch viscosity model. 

2 - Viscosity is Input 

LW 

tl 

ft 

Specifies the laminar viscosity model used for the 
wall cooling modes. LW = 1 employs the 
Sutherland air viscosity model. 
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O /H 2 characterization program input parameters 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

NTYPE 

Combustion 

I 5 

0 - Axisymmetric flow field 

1 - Plane two-dimensional flow 

ISOBAT 

ft 

tt 

1 - Isoenergetic wall 

2 - Regene ratively cooled wall 

3 - Transpiration cooling 

4 - Slot cooling 

MB 

ft 

ft 

0 - Free jet; P (x) prescribed 

1 - Ducted flow; P (x) prescribed 

2 - Ducted flow; Wall radius (x) prescribed 

ICHEM 

ft 

tt 

1 - Chemically frozen flow 

2 - Equilibrium ("complete combustion") 

chemistry 

3 - Finite rate chemical kinetics 

MC 

tt 

tt 

Printout of the flow field radial profiles is made 
every MC finite difference steps . The default 
is 10. 

MG 

f T 

ft 

Specifies the diluent specie used as a tracer in the 
transpiration model. The diluent may also be 
present initially in the main stream flows, when 
other wall cooling models are used. 




1 - Diluent is Nitrogen 

2 - Diluent is Helium 

3 - Diluent is Argon 

LZ 

MA 

MY 

MH 

tt 

tt 


Printout dump controls for various portions of the 
program. In general, nn = 0 means no dump; 
nn = 1 yields moderate dump; nn ^ 2 yields 
overwhelming, printout dump. . 

ISBATY 

tf 

tt 

Specifies alternatives in the regeneratively 
cooled and transpiration cooled wall models: 




Regenerative Cooling (ISOBAT = 2) 

ISBATY 

1. Wall temperature , T w (x) , (°R), specified. 

The wall heat transfer, q (x) Btu /in 2 sec), 
is computed. 


q w (x) specified. T w is computed. 


2 . 


S-1220 


o /h characterization program input parameters 

(Continued) . 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

ISBATY 
(Cont. ) 

Combustion 

I 5 

Transpiration Cooline (ISOBAT = 31 
ISBATY 

1. Coolant temperature, T c (x), (°R), and coolant 

unit area mass flow rate, (Ou (x) ) (lb/in^sec), 

c 

are specified, and coolant domain "edge" con- 
ditions are used in the model. T w is com- 
puted. 

2. T c (x) and (Ou (x) ) c are given, and bulk flow 
"edge" conditions are used. T w is computed. 

3. T c (x) and T w (x) are given, with coolant domain 
"edge" conditions. ( Pu) is computed. 

4. T c (x) and T w (x) are specfied with bulk flow 
"edge" conditions. (Ou) is computed. 

NS LOT 

ft • 

It 

The number of slots in the slot wall cooling model. 

PRNT 

It 

E10.8 

Printout interval Ax (inches). 

XMAX 

" 

M 

Axial distance to which calculation is carried out, 
Xmax ( in ches). 

X 

M 

" 

Axial station at which calculation is begun, 
x initial (inches) . 

XLE 

It 

tf 

Turbulent Lewis number. 

SIGMA 

II 

IT 

Turbulent Prandtl number. 

DELPSI 

If 

• tf 

The spacing between adjacent flow field grid points. 

The dimensions of A # are / lb/sec in axisym- 

metric flows and (lb/sec) in plane two-dimensional 
flows . 

XMPS 

1 1 

ft 

Available for use in reducing the diffusion step 
size, Ax = Ax /XMPS. Let MPS = 1. 

XK2 

tt 

TT 

Constant employed in the turbulent viscosity models. 

PSI(l) 

II 

ft 

Value of the flow field grid point, ip 1, nearest the 
chamber centerline , with dimensions of 

£ lb/sec j( ™PE). 
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oVH 2 characterizationprogram input parameters 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

XK(1) 

Combustion 

E10.8 

ol in the Hersch viscosity model. 
Let CL- 1. 

XK(3) 

M 

ft 

The spacing of initial fuel and oxidizer "rings" in 
the Hersch model, s (inches). 

XK(4) 

tf 

ft 

Parameter in the Hatch-Papell film cooling model. 
Input as . 04. 

XK(5) 

tf 

f f 

Mean maximum number of finite rate chemistry 
steps per flow field diffusion step. Default value 
is 10. 

XK(6) 

t» 

tf 

C f /2 in the momentum eq. wall boundary con- 
dition. Recommend using 1x10 . 

XK(7) 

M 

tf 

D, for x/D printout, (inches). 

TAR 

ft 

ft 

Initial wall radius (inches) for MB = 0 and 1; 
initial static pressure (lb/in ) for MB =2. 

XP(1) thru 
XP(4) 

It 

ft 

End points of domains of polynomials for wall 
radius, or static pressure, and other wall 
boundary conditions (inches) , which are im- 
put below. 

XP(5) 

ft 

tf 

Maximum finite rate chemistry time step(seconds) . 
Default is 1 x 10"^ seconds. 

XP(6) 

If 

ft 

Lower tolerance for changes in finite rate chemis- 
try time step. Default is 5 x 10 

CGP(I, J) 

ft 

7E10.8 

Coefficients of polynomial J for static pressure 
(MB=0,1) or chamber wall radius (MB=2). 

J =1,4 (four cards) . 

6 '-I 

F j( X) = i f i a j ( X ' X I*) J for x < xp d) 




where CGP (I,J) = X* 

TWX(I, J) 

ft 

ft 

Four sets of polynomial coefficients for T w (X) (°R) 
For regenerative cooling if ISOBAT = 2. 

For transpiration cooling if ISOBAT = 3. 

QLX(I, J) 

f f 

tt 

Four sets of polynomial coefficients for q w (X) 
(Btu/ in 2 gec ) (Regenerative cooling) . 


S-1220 


o £ /h 2 characterization program input parameters 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 


RUCX(I,J) 

Combustion 

7E10.8 

Four sets of polynomial coefficients for ( pu) c 
(lb/in 2 sec) (Transpiration cooling). 

TCX(I, J) 

tf 

If 

Four sets of polynomial coefficients for T c (X) (°R) 
(Transpiration cooling). 

NS LOT 

ft 

E10. 8 

Number of slots in the chamber wall ( ^ 21 ). 

XS(K) 

ft 

11 

Axial location of first slot (inches). \ 


SH(K) 

t! 

11 

Height of first slot (inches) 

1 

UC(K) 

11 

11 

Coolant velocity (in/sec) 

>K = 1,2 

RUCF(K) 

11 

11 

2 1 

Coolant mass flux (lb/in sec) | 

NSLOT 

TCS(K) 

11 

11 

Coolant temperature (°R) > 


* RSTAR 

11 

11 

Radius of curvature of engine throat (inches) 

PST 

11 

11 

Throat diameter (inches) 


PRC 

M 

11 

Coolant Prandtl number. 


NELEM 

Injection 

I 3 

No. of injector elements 


DCHAMB 

11 

F10.0 

Chamber diameter at injector plane (inches) 

TO 

11 

ft 

Oxygen total temperature (°R) 


TH 

11 

It 

Hydrogen total temperature (°R) 


EMR 

11 

11 

Rupe mixing factor 


WT 

11 

It 

Mass flow through injector (lb/sec) 


PO 

11 

If 

Chamber total pressure (psia) 


FM 

11 

ft 

Mass fraction oxygen thru injector 


T(I) 

Combustion 

E10.8 

Static temperature (°R) at grid point 

I 

U(D 

T? 

M 

Axial velocity (in/sec) at each grid point 

FIX(I) 

H 

ft 

Chemistry time step (seconds) at each grid point 

TELA P (I) 

11 

ft 

Integrated streamline residence times (seconds) 
for each grid point 

YSPEC(l.I) 

11 

f 1 

H mass fraction at grid point I 


YSPEC(2,I) 

11 

If 

O mass fraction at grid point I 


YSPEC(3,I) 

11 

ft 

HgO mass fraction at grid point I 
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O /H,. characterization program input parameters 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

YSPEC(4 , 1) 

Combustion 

E10.8 

H mass fraction at grid point I 

YSPEC(5 , 1) 

IT 

ft 

O mass fraction at grid point I 

YSPEC(6 , 1) 

If 

If 

OH mass fraction at grid point I 

YSPEC(7 , 1) 

If 

It 

HO mass fraction at grid point I 

YSPEC(8 , 1) 

ft 

II 

H O mass fraction at grid point I 

Z Z 

YSPEC(9 , 1) 

IT 

If 

Diluent mass fraction at grid point I 

RT 

Performance 

F10.0 

Throat radius (inches) 

RE 

II 

f f 

Exit radius (inches) 

XN 

If 

If 

Nozzle length (inches) 

PERBEL 

ft 

tf 

Percent bell 

ENTHO 

IT 

If 

Oxygen enthalpy (Btu/lb) 

ENTHH 

ft 

t r 

Hydrogen enthalpy (Btu/lb) 

OFINPT 

IT 

ft 

Mixture ratio 

PER 

ft 

ti 

Percent fuel injected in supersonic region 

TTT 

It 

ft 

Temperature of fuel injected in supersonic region 

NNODE 

Heat Transfer 

I 3 

Number of nodes 

NTYPFL 

1 1 

If 

1 if transient, 0 if steady state 

NFLMFL 

ft 

It 

1 if film cooling, 0 if not 

NRGNFL 

i r 

II 

1 if regen cooling, 0 if not 

NLFL 

i r 

II 

1 if liner, 0 if not 

INJFL 

IT 

ft 

1 if injector, 0 if not 

TSTOP 

I! 

F12.0 

Cut-off time for transient case (sec) 

NCOUNT 

II 

If 

Number of iterations for steady state case 

PRINT 

INTERVAL 

ft 

tl 

(input as a real number) 

Print interval in seconds for transient and in number 

TWA LLF 

ft 

f 1 

if iterations for steady state 
Initial wall temperature (°F) 


S-1220 


°2^ H 2 CHARACTERIZATION PROGRAM INPUT PARAMETERS 

(Continued) 


NAME 

MODEL 

FORMAT 

DESCRIPTION 

TSINKF 

Heat Transfer 

F12.0 

Sink temperature (°F) 

CAPPA 

tt 


Conductivity (Btu/in sec °R) 

EPS 

ft 

TT 

Wall emissivity (none) 

PO 

ft 

TT 

2 

Combustion chamber pressure (lb/in ) 

OF 

Tt 


Mixture ratio (none) 

RC 

tt 


Radius of curvature at throat (inches) 

RHO 

n 

TT 

3 

Density of wall material (lbs/in ) 

CP 

f T 

f T 

Specific heat capacity of wall material 
(Btu/lb °R) 

TINJH2 

TT 

M 

Hydrogen injection temperature (°R ) 

TIN J02 


- 

Oxygen injection temperature (°R) 

X 

tt 


Axial distance of node (in) 

DI 


ft 

Wall inner diameter (in) 

DO 

TT 

TT 

Wall outer diameter (in) 

NINJ 


I 3 

Number of film injection stations 

SLOT 

TT 

F12.0 

Axial distance of film injection station (in) 

HSLOT 

TT 

ft 

Slot height (in) 

WCOOL 


IT 

Coolant weight flow (lb/sec) 

HRWD 

T f 

ft 

Hydrogen regen weight flow (lb/sec) 

NPASS 

f f 

I 12 

Number of regen passages 

XREGEN 

Tt 

F12.0 

Regen injection station (in) 

X 

TT 

TT 

Axial distance of node (in) (must correspond to 
nodes) 

HPASS 

ft 

ft 

Regen passage height (in) 

WPASS 

Tt 

TT 

Regen passage width (in) 

FIN EFFICIENC 

Y 

ft 

Fin efficiency (none) 

MASS FLOW 

Tt 

ft 

Mass flow between liner and wall (lb/sec) 

EMISSIVITY 

ft 

rt 

Liner emissivity (none) 

LINER LENGTH 

ft 

M 

Maximum axial distance of liner (in) 
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O z /H 2 characterization program input parameters 

(Continued) 


NAME 


MODEL 


FORMAT 


DESCRIPTION 


X 

LINER ID 

LINER OD 

AINJS 

EMINJ 

ARINJH 

HGINJH 

ARINJO 

HGINJO 

RESINJ 

OFINJ 

WTINJ 

CPINJ 


Heat Transfer 


F12. 0 


Axial distance of node (in) (must correspond 
to a node) 

Liner inner diameter (in) 


Liner outer diameter (in) 

2 

Injector surface area (in ) 

Injector emissivity 

2 

Injector H Convection Area (in ) 

Li 

Injector H 2 Convective Coefficient 


H *- ) 

\ m sec °R/ 


Injector O Convection Area (in ) 

Li 


Btu 


Injector C> Convective Coefficient/ T~i 

2 yin z sec “R, 

Injector - Combustor thermal resistance 
(sec°R/Btu) 

2 

Injector face heat flux (But/in sec) 

Injector weight (lbs) 

Injector specific heat (But/lb °R) 
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SAMPLE CASES 
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S-1220 


Sample Case I 


Sample Case I is a combustion and performance computation utilizing the injection 
model to provide starting profiles. The first data card required is the run control card 
with NCOMBU and NINJE set equal to 1 (all others zero or blank). All dynamics input 
is omitted. The second data card, then, is the combustion program title card followed 
by the combustion control card (MPSI thru NSLOT). The fourth, fifth, and sixth cards 
are (PRNT thru X), (XLE thru PSI(l) ) , and (XK(1) thru XK(7)). The next five cards 
define the wall contour. The seventh card (TAR thru XP(6) ) defines the polynomial limits, 
and the next four are the polynomial coefficients (CGP(1, J) thru CGP(7, J) for J= 1, . . ,4). 
The next data cards are the two injection model cards, (NELEM thru EMR) and (WT thru 
FM). All intervening cards shown on the input sheet are omitted. The final cards are the 
two performance input cards (RT thru ENTHH) and (OFINPT thru TTT). 

The following page is a copy of the Sample Case I input report printed by the 
program followed by a sample of the combustion output and the final performance report. 
Note that injection model results (species, velocity, and temperature profiles) are printed 
as part of the combustion input report. 


100 


TEST Case AEROJET H2/02 1500 U3 thrust engine 
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the mixture ratio is a.ib 

The I SR AcTeR COMEUSTiON BUT WIThCUT NOZ2LE LOSSES IS • «07,&5 
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Sample Case II 


Sample Case n is a dynamics computation for the baseline Aerojet engine. The 
first data card required is the run control card with NDYNA = 1 and NBLFG = 1, all 
others zero. This card is followed by the dynamics control card (NNODE thru CD). 
Since a blowdown calculation is specified by NBLFG = 1 , the next card defines the blow- 
down volumes (NBLOW(I)). Then, NNODE cards are required to define the volumes 
followed by NCONN cards to define the flow passages. The last two cards are the 
ignition cards (SPTIME thru SPKF) and (SPGAP thru DC). 

On the next page is a copy of the Sample Case n input report followed by several 
samples of the dynamics output. 
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• • • DYNAMICS INPUT • • • 


NNOCE 

NCONN 

NITER 

nchek 

PCHEK 





7 

8 

4000 

1 

1,000 





NVOL 

I COME 

VOL 

PRES 

TEMP 

RM1X 




1 

1 

48.300 

.000 

530.000 

,000 




2 

-o* •* 


.000 

530.000 

.000 




3 

-o 

23.800 

.000 

530.000 

.000 




4 

-o 

32.800 

.000 

53C.00U 

.000 




s 

3 

. 4 A 7 

.000 

530-000 

.000 




6 

'CllQODOO-OUO 

375.000 

530.000 

.000 




/ 

-oioo 

000.000 

375.000 

530.000 

1,000 




I ADMIT 

I ADM I T 

! RTvpjr 

CARE A 

CCOEFP 

TJMON 

TIMOFF 

T I MQPN 

T I MCL.S 

1 

2 

2 

2.395 

1.000 

-.0000 

*.0000 

-.0000 

-.0000 

5 

1 

0 

.035 

l-ouQ 

-.0000 

-.0000 

-.0000 

-.0000 

3 

1 

0 

,438 

1.000 

-.0000 

-.3000 

*,0000 

-.0000 

4 

1 

0 

,34V 

1.000 

-.0000 

-. 3000 

-.0000 

-.0000 

6 

3 

0 

.500 

1.000 

. 0280 

. 1300 

.0100 

.0087 

7 

4 

0 

,332 

1.000 

.0260 

. 1223 

.0087 

.0074 

6 

5 

0 

,0U7 

1.000 

.0100 

. 1300 

.0100 

.0100 

7 

5 

0 

.012 

l.UOO 

.ocou 

, 1200 

.0100 

,0100 


SPTJ*E SPrfT L 
, 0 U 0 C ,C50o 


SPKP 

,0u25 


SPGAP 

,u 2 Sc 


SP4RKP SPaRkE 
5Q0L o . o lOO.u'OCD 


DC 

3. 6Q00 



•• DYNAMICS RESULTS *• 
TIME ■ ,026193 


ITERATION NO 6C0 




VOL NC 

PRESSURE 

TEMP 

( cx 


< p Si A ) 

(R > 


1 (CCN8USTCP) 

13-82 

5226,4 

93 . Bpi 

2 

• Cl 

53«0.7 

89,020 

3 

13.55 

4146,1 

88,917 

4 

42.46 

1557.7 

99.031 

5 (PILOT) 

295 . V6 

6376.4 

86.754 

6 

374.99 

530.0 

,000 

7 

374.99 

530,0 

100.030 

DEiT L5E0 -000025 



TuPLST* 

,00 

IMPULSE* 

,00000 

ACCUM V-T FLOW* 

• CiU-'CQ 



VOL 1 PPES 

CHANCE* 

4.93GPSI 

PER MS 
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•• DYNAMICS RESULTS •• 
TIME * ,028671 


ITERATION NO 
VOL NC 


620 


(CCmBUSTOR > 


1 
2 

3 

4 

5 

6 
7 

oelt lsed 


(PILOT) 


•000025 


THRUST* 

ACCUM »T FLOV>» 


RESSUPE 

temp 

t ox 

(PSlA) 

< R ) 


17.16 

5190.5 

94,218 

•Cl 

5368,6 

89.456 

18.43 

1531.5 

61.420 

5J.91 

1351,5 

99,245 

295.96 

6376,4 

86.754 

374.99 

530.0 

,000 

374.99 

1 

530,0 

100.000 

)0 

impulse* 

, 0000 c 


VOL 


.oocco 

pR£S CHANGE* 


8.043PS1 PER MS 


ISP* 


,00 


•• DYNAMICS RESULTS •• 
TIME * ,029144 


ITERATION NC 640 




VOL NC 

PRESSURE 

temp 

( OX 


(PSIAI 

IR) 


1 (CCMSUSTOP) 

22 • C8 

5434,0 

92,102 

2 

• Cl 

5363.6 

89,829 

3 

26.41 

697,3 

35.323 

4 

59. «6 

1199,5 

99.396 

5 (PIUCTJ 

295.95 

6376,4 

06. 754 

6 

374.99 

530.0 

.000 

7 

374.99 

530.0 

100. oco 

DECT USED .000025 



THRUST* 

,CU 

IMPULSE* 

.00000 

ACCUM WT FLO'** 

.uooco 



VOL 1 PFES 

CHANGE* 

11.550PSI 

PER MS 


•• DYNAMICS RESULTS •• 
TIME » ,029656 

:TERATION NO 66q 


VCL NC 

PRESSURE 

temp 

( ox 


< PS I A ) 

IR> 


1 (CCMBUSTOR) 

29.<4 

5579,8 

88.330 

2 

•Cl 

5380. 3 

89,878 

3 

36.97 

701,9 

20,038 

4 

66.48 

1084,6 

99.505 

5 (PILCT) 

295-95 

6376,4 

86.754 

6 

374.98 

530,0 

.000 

7 

374.99 

530.0 

100,000 

OElT LSEO .000025 



THRUST* 

,0C 

IMPULSE* 

.00000 

ACCUM N'T Fl.04* 

.OUC'CO 



VOL 1 P p £S 

CHANCE* 

26 .78VPSI 

PER MS 
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Sample Case HI 


Sample Case HI is a steady state heat transfer computation for the baseline 
Aerojet engine using a combination of film and regenerative cooling. The first data card 
is the run control card with NHEAT = 1. Next comes the heat transfer title card followed 
by the two control cards (NNODE thru INJFL) and (TSTOP thru PRINT INTERVAL). 

The next two cards are (TWALLF thru OF) and (RC thru TINJ02). The next NNODE 
cards define the combustor geometry. Then come the film cooling cards; the first 
specifies the number of injection points and the following NINJ ones define injection 
locations, slot heights, and coolant mass flows. The last data cards required to run 
Sample Case M define the regenerative cooling. One card (HRWD thru XREGEN) 
specifies the regenerative coolant mass flow, the coolant introduction point, and the 
number of coolant passages. The remaining cards specify coolant passage geometry at 
each nodal point. 

On the next page is a copy of the Sample Case El input report followed by several 
samples of the heat transfer output. 
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—ROCK ET 

HF&T TP ANSWER 

MODEL- 

AEROJET 

REGENERAT1VELY 

CCCLEO C 


ITERATION NO. CORE 

0/F 

CORE TEMP 

WT. F LCW 



? 




DEC F 

LR/SEC 





70 4. 

OCO 

5185. 3 

3.057 



1 

NOOE 

STATION T 

WALL 

T FILM 

T BULK T 

LINER 

T H2 

A 

> 



- 



REGEN 


LINER 

1 

* 

1 

0 .0 

141.5 

141.4 

141 .4 

0.0 

0 .0 

It 

2 

0.01 

141.5 

144.3 

141.4 

0.0 

0.0 


3 

0.65 

144.3 

308.2 

141.4 

0.0 

0 .0 

• I 

4 

1.30 

147. 1 

472.5 

140. 8 

0.0 

0.0 

:i 

5 

1.95 

149.5 

633.9 

135.6 

0.0 

0.0 

if 

it 

6 

2.60 

151.6 

751.0 

137.7 

0 .0 

o.o 


7 

3.25 

155.4 

546.6 

135.5 

o.c 

0 .0 

-i 

e 

3.90 

162.3 

1111.5 

132.8 

0.0 

0 .0 

it 

9 

4.55 

174.9 

1207.5 

129.5 

0.0 

0.0 

it 

jt 

to 

5.20 

193.2 

15C1.3 

125.2 

0 . 0 

0.0 

il 

M 

u 

5.80 

158. 0 

1699.7 

115.6 

0.0 

0.0 

II 

12 

6.88 

182.2 

2016. 9 

105.9 

0.0 

0.0 

11 

12 

7. <16 

14Q.0 

2 2 ''o . 2 

37.7 

0.0 

0.0 


14 

9.04 

168.0 

2461.5 

70.0 

0.0 

0 .0 

V 

n 

15 

9.05 

168. 2 

143.3 

0.0 

0.0 

0.0 

M 

16 

10.12 

216. 1 

327.7 

0.0 

0.0 

0.0 

11 

17 

11.20 

337.5 

500.3 

0.0 

0.0 

0.0 

)l 

18 

12.27 

284.5 

661. 7 

0.0 

0 .0 

0.0 

11 

19 

12.28 

284.5 

656.2 

0.0 

0.0 

0.0 


— ROCK FT HEAT TRANSFER MODEL — AFPCJET PE GE NE P AT ! VE IT CCCIEC C2/H2 ENGINE 



ITEPAT ION 

NO. CORE 0/F 

COPE TEMP 

WT. FLOW 



I 




OEG F 

LB/SEC 





50 

4.000 

5185. 2 

3.057 



, ’ 

NODE 

STAT ION 

T WALL 

T F I 1.m 

T BULK T 

LINER 

T H2 

» 





REOFN 


LINER 


1 

O.C 

140.5 

140.8 

140.9 

0.0 

0.0 


2 

0.01 

140.5 

143.5 

140.9 

0.0 

0.0 


3 

0.65 

143.8 

3C 7. E 

140.9 

O.G 

0 .0 

/ 

4 

1.30 

146.7 

4 7 2.6 

140.4 

0,0 

0.0 


5 

1.95 

149.1 

633.6 

13 c .2 

0.0 

C .0 


6 

2.60 

151.2 

750.8 

137.3 

O.C 

0.0 


7 

3.25 

1 55.0 

C46.5 

135.1 

o.o 

0.0 


8 

3 . °0 

161.9 

1 110.4 

132.4 

0.0 

0 .0 


Q 

4.55 

174.5 

12=7.5 

125.0 

0 .0 

0.0 

IK 

10 

5.20 

152.8 

15 Cl. 1 

124.8 

0 .0 

0 .0 


u 

5.eo 

157.6 

1690.9 

119.2 

0.0 

C .0 


12 

6. se 

iei.8 

2017.2 

105.5 

0.0 

0 .0 


13 

7.96 

147. 1 

2274.6 

87.3 

o.o 

0.0 


14 

9.04 

146.5 

2462.3 

70.0 

0.0 

0 .0 


15 

9.05 

146.7 

142.7 

0.0 

. 0.0 

0.0 

19 

16 

10.12 

158.3 

327.3 

C.O 

0.0 

0 .0 


17 

11.20 

316.4 

500.0 

0.0 

0.0 

0.0 


13 

12.27 

237.8 

661,5 

0.0 

0.0 

0.0 


19 

12.28 

237.8 

656.0 

0.0 

0.0 

0.0 
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—ROCKET FEAT TRANSFER Mf.OEl— AEROJET R EGENE R AT I VE l Y CCCLEfr G2/H2 ENGINE 


ITERATION NO. 

CORE 0/F 

CORE TEMP 
OEG l 

WT. nCV 

LO/SEC 


60 

A .000 

5185.3 

3.057 


NOOE STATION 

T WALL 

T F 1 LM 

T 31JLK T LINER 
RFOEN 

T H2 
L INER 


i 

n.o 

1 A l . 1 

1A1.1 

1 A 1 . 2 

0.0 

0.0 

2 

0.01 

1*1.1 

l *3. a 

I A 1 .2 

0 .0 

C .0 

3 

0.65 

1 A* . 1 

30R.0 

1 A 1 . 2 

0.0 

0.0 

4 

1 . 30 

1 46. 9 

472. 8 

1A0.6 

0.0 

0.0 

5 

1.05 

1A5. 3 

633. 8 

139. A 

0.0 

0 .0 


6 

2.60 

151. A 

750.9 

137.5 

0.0 

0.0 

7 

3.25 

155.2 

<546. 6 

135.3 

0.0 

0.0 

8 

3.90 

162. 1 

1110.5 

132.6 

0.0 

0 .0 

9 

A. 55 

1 7A. 7 

1297.5 

129.3 

0 .0 

0.0 

10 

5.20 

193.0 

1500.0 

1 25.0 

0.0 

0.0 

“ 11 

5.80 

157.8 

1695.8 

115. A 

0.0 

0.0 

12 

6 . P 8 

182.0 

2C1 7.1 

105.7 

0.0 

0.0 

13 

7.96 

1A8. 1 

22 7A.A 

87.5 

o.c 

0.0 

1A 

9.0A 

157.8 

2A62.1 

70.0 

0.0 

0.0 

15 

9.05 

158.0 

1 A 3 . 0 

0.0 

O.P 

0.0 

“ 16 

10.12 

207.7 " 

327.6 

0.0 

o.c 

0.0 

17 

11.20 

327. A 

5C0.2 

0.0 

0.0 

0.0 

18 

12.27 

262.3 

661.6 

0.0 

0.0 

0.0 


12.28 

262.3 

656.1 

0.0 

0.0 

0.0 

.. 

... 

- 

.. 
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SUBROUTINE DESCRIPTIONS 


MAIN 


The MAIN routine accepts input for all models except heat transfer. Unit 
conversions are made in MAIN to prepare input data for the combustion subprogram 
which makes calculations in metric units. The MAIN routine also contains the logic 
for calling the principal models. It controls the dynamics model shutdown, time ad- 
vance, etc. 

SUBROUTINE AREA 

When the thrust chamber wall radius, r w (x), is specified (MB = 2), dP/dx is 
evaluated in this subroutine by means of an iterative process involving the strict conver- 
gence of r 7 w = r' w (dP/dx, physical chamber conditions at x) to r w / (x). 

SUBROUTINE BEGIN 

Molecular weights, heats of formation, and polynomial curve fits for the enthalpy 
and entropy of the chemical species are stored here. In addition, the units of the thermo 
dynamic property information are converted from the English and c.g.s. systems to the 
S. I. system. 

SUBROUTINE BLOW 


This subroutine is used to calculate chamber pressure after the valves are closed 
Its purpose is to reduce the computer time required to calculate blowdown conditions. 

SUBROUTINE BUMP 


BUMP is used to reduce the dynamics calculating increment during the first 
few iterations when the stable time is relatively large. The purpose is to increase 
the resolution of the model during the first iterations. 

SUBROUTINE BULK 

- i ^ lr • . . 

This subroutine integrates the dependent variable (u, H, and a J or at ) 
radially across the chamber, at each axial station, to determine their bulk values, and 
then calls the appropriate subroutines to determine bulk values for the dependent vari- 
ables (h, T, Cp, p, w, etc.). These bulk quantities are used for printout purposes , 
and where needed, in applying the chamber wall boundary conditions. 
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SUBROUTINE CFLOW 

CFLOW is used to calculate mass flow for a choked duct given the 4fL/D 
parameter. 

SUBROUTINE COLD 


Subroutine COLD calculates performance of engine which propellants are not 

ignited. 

SUBROUTINE COMB 

COMB is the calling program for the combustion module. Once called by MAIN, 
COMB controls combustion calculations much as MAIN controls the overall characteri- 
zation run. 

SUBROUTINE CONSRV 

An explicit solution of the conservation equations for velocity, total enthalpy, and 
either the chemical specie mass fraction (frozen (ICHEM = 1) or finite-rate chemistry 
(ICHEM =3), or the chemical element mass fraction (equilibrium (ICHEM = 2) chemistry) 
diffusion equations is performed in this subroutine. 

SUBROUTINE CPRF 

This routine calculates the choking pressure ratio for a duct for a given 4fL/D 
parameter. 

SUBROUTINE CRSPLT 

Crossplot data from GRAD against composite kinetic curve of NASA CR 72601 
to find Freeze Area Ratio. 

SUBROUTINE DENSE 


The density and molecular weight at each grid point in the flow field are computed 
here as a function of p(x), T (x, ip) and the cd (x, ip ). 

SUBROUTINE DIV . BL. SIDEL 

DIV calculates percent I S p loss due to divergence. 

SUBROUTINE BL 

Subroutine BL calculates boundary layer loss. 
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SUBROUTINE SIDEL 

SIDEL converts percent recombination to specific impulse. 

SUBROUTINE DYNAM 

DYNAM contains the logic for calculating mass flow rates and volume mixtures, 
temperatures, and pressures. It contains the logic for sparking, ignition, and 
quenching tests, and also calculates the stable time and the instantaneous orifice cross- 
sectional areas . All dynamics output information is printed from DYNAM. 

SUBROUTINE EQUILC 

When the option of equilibrium chemistry is chosen (ICHEM = 2), the chemical 
species composition at each flow field grid point is computed here using a "complete 
combustion" model. The specie mass fractions ( a- 1 (x, 0) ) for the species H 2 0, H 2 
and O 2 are computed from the element mass fractions (3(x, 0)) for H and O. 

SUBROUTINE FFCF 

Subroutine FFCF is used to compute duct friction coefficient given the flow 
Reynold's number assuming a smooth duct surface. 

SUBROUTINE FG 

The coefficients and forcing vector of the linearized chemical kinetic reactions 
are evaluated in this subroutine. 

SUBROUTINE FGX 

The coefficients and forcing vector of the linearized chemical kinetic reactions 
are evaluated in this subroutine. 

SUBROUTINE FLOW 

FLOW is used to calculate mass flow rate in ducts where choking does not occur. 
Two parameters are required, the pressure ratio across the duct and the friction 
parameter, 4fL/D. 

SUBROUTINE FLUX 

When either the transpiration-cooling model, or film-cooling model, is being used 
at the wall, this subroutine monitors the amount of mass being added to the flow field, and 
adds additional grid points to the computation as required. 
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SUBROUTINE FMOLWT 

This routine calculates the molecular weight of the products of GO^/GH^ 
combustion given the mixture ratio. 

SUBROUTINE GRAD, NOZ 

GRAD determines the area ratio gradient of the nozzle at specified area ratios 
based on NASA CR 72601. NOZ is used to calculate dA/dX from NASA CR 72601. 

SUBROUTINE GRID 

This subroutine controls the addition, or subtraction, of streamline grid points 
from the finite difference flow field computation. Grid points are added above (until 
the wall is reached) and below (for ^(1)> 0) the present grid according to tests involving 
the principal flow field variables (H, u, and oP or 3 ). In addition, when one less than 
double the initial number of grid points is in use, alternate grid points are discarded, the 
interval between grid points doubled, and the computation returns to using the initial 
number of grid points. 

SUBROUTINE HEAT 


The thermodynamic properties at each flow field grid point are computed here. 
There are two options. For KOPT =1, C p (x, ip), h (x, <p), and H (x, ip) are determined 
from T (x, ip), u (x, ip) and Cf> (x, ip). For KOPT = 2, T (x, ip), C (x, >P), and h (x, ip) are 
determined for H (x, ip), u (x, ip) , and Of 1 (x, ip). P 

SUBROUTINE HE ATT 

HEATT is the main heat transfer routine. It accepts heat transfer input and 
creates the thermal model. It contains the logic to calculate temperatures using any of 
the available cooling models. HEATT also prints results and input parameters in a 
concise format. 

SUBROUTINE HONC 


The finite rate chemistry calculation is controlled in this subroutine. It regulates 
the chemistry time steps, calls the subroutines which (a) compute the reaction rates, 

(b) evaluates the linearized equations using a Pade' rational approximation for the ex- 
ponentials representing the solution of the coupled first order linear ordinary differential 
equations , and (c) solve the matrix which represents the integration of the chemical 
kinetics equations over the particular time step. 

SUBROUTINE INJECT 

Subroutine INJECT is used to compute starting profiles for the combustion model. 
It calculates velocity, temperature, and species mass fractions at each combustion grid 
point. It zeros all species mass fractions other than those of 0 2 and H 2 . 
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SUBROUTINE INVERT 

This subroutine serves the same purpose as HEAT, but for a single grid point 
at a time rather than the entire grid. "Entry HOOT" is part of this subroutine. 

SUBROUTINE MARCH 


This subroutine performs the calling function for the combustion calculations. 

It calls the subprograms to make sequential calculations in the axial direction. 

SUBROUTINE PAPE 

This subroutine controls the matrix solution of the integrated linearized chemical 
kinetics equations and prepares the inputs for the actual matrix calculation. 

SUBROUTINE PERF 


This subprogram is the main driver for computing the performance of the rocket 
engine using data transferred from the combustion and dynamic programs. Igp, Thrust, and 
Cp are calculated. 

SUBROUTINE PERSHF 

This routine determines percent recombination as a function of Freeze Area Ratio 
from CRSPLT. 

SUBROUTINE PREPAR 


When the finite rate chemistry option (ICHEM = 3) is employed, this subroutine 
serves as a connecting link between the flow field computation and the chemical kinetics 
model. Here inputs are prepared for the chemistry subroutines, the chemistry time steps 
related to the flow field step size, and the results of the chemistry computation prepared 
for insertion into the flow field computation. 

SUBROUTINE PRESS 

This subroutine is used to evaluate an arbitrary function and its first derivative, 
f(x) and df (x)/dx, from input polynomials of as high as fifth order. 

SUBROUTINE PICRT ' ' 


The solution of the matrix representing the linearized chemical kinetic equations 
is performed here. 
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SUBROUTINE PRINT 


This subroutine is used to printout the radial profiles of the principal flow variables , 
such as H, u, t, a? , h, etc. In addition, dependent variables which are not required for 
the flow field computation, but are of interest, such as Cp, 7, <P, H s , etc. , are computed 
and printed. Also, average properties across the flow field are computed and printed. 

SUBROUTINE RTNI 

RTNI is a routine which used Newton's method to solve nonlinear equations. It 
is used by INJECT to solve for flow Mach number. 

SUBROUTINE SI 

Subroutine SI calculates equilibrium and frozen specific impulse for each 
streamline. 

SUBROUTINE SLOT 

When slot cooling is used at the chamber wall (ISOBAT = 4), this subroutine is 
used to apply the wall boundary conditions to the conservation equations for energy and 
diffusion. 

SUBROUTINE SPHEAT 

SPHEAT calculates the ratio of specific heats of the products of G0 2 /GH 2 
combustion given the mixture ratio. 


SUBROUTINE STEP 

The step-size, Ax, for the next flow step is computed here. The step-size is 
defined by applying the von Neuman stability criterion to the conservation equations. This 
results in the criterion that 


Ax = minimum 



(A !)) 2 \ 

4jiSc / 


(£$) 


m,n' 1 \ 

(— ) / 

JL\ 

\Sc‘n,Ta+h ' 

Scl 


' -r ) 
,m+l n,m 


SUBROUTINE TGR 

TGR is used to calculate equilibrium combustion gas temperature for given mixture 
ratio and pressure. It is utilized in the dynamics model and the heat transfer model. 
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SUBROUTINE TRANSP 

When transpiration cooling is used at the chamber wall (ISOBAT = 3), this subrou- 
tine is employed to apply the wall boundary conditions to the conservation of energy and 
diffusion equations . 

SUBROUTINE VASC 

This subroutine calculates a laminar viscosity, using the Sutherland model for 
air, for use in applying the wall boundary conditions. 

SUBROUTINE VISC 


The turbulent viscosity coefficient, p. (x), is computed here using a model employ- 
ing Hersch's mixing parameter: 

ak Dpu 
K 1+x/S 


SUBROUTINE WALL 

This subroutine controls the application of the wall boundary conditions to the flow 
field computation. After first testing as to whether the wall lies above, or is in the flow 
field at a given x station, it proceeds to apply the boundary conditions to the momentum 
conservation equation, and then calls the appropriate subroutine to apply the B. C. to the 
energy and diffusion conservation equations. 
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SUBROUTINE FLOW CHARTS 
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SUBROUTINE AREA 


/this the\ 
'first axial ' 1 
step in the 
^romputatioiy 
\ ? / 


NUMBER - 0 
LT 40 
LU = 20 
LX ^ 20 


Initialize parameters for the 
computation 

\dx/n + 1 from 
the chamber mass flow, area, 
and flow field properties. 

— ► 





NO / Is \ 
— < PNEW > . 0 N 


Use first guess for 


Vdx / n 4 . i 'dx/ n 
PNF.W - P n 


Call 'PRESS” to evaluate 
the wall radius, y (x n v j ) , 

from the inputted polynomials 



(DO Loop) 

1 

Compute Rr V 


PNEW = 

G7 x P n 

(dfi.S = 

PNEW - P n 


dx 

HDPDX = 



'dx/n+ 1 



' WRITE 

YW.DPDX.PNEW, 
P.DX.X / 


DO LOOP 
M = 1 to LT 


• PNEW is too large; readjust the iteration procedure 


IOl'T, I, LR, MPSI.KAT 


PNEW, RU(I), RUR(I), DX 


DPDX, RUT(I) 


PFRD = PNEW 


&U 1 

UFAC = 

PNEW - 
dx 

dxx (^) 1 
»*i 


l?u ( ’l 


PNEW = . 9 *PFRD 



is still large 
-.enough to 
cause numerical 
errors in momentum 
computation - readjust: 

PFAC “.lx u'ejiu ) 1 

_ Q. 

dx 


PFAC = - PFAC 


* PFAC 

dx 'n - 1 

PNEW = P + dx x 
a 

KAT = KAT + 1 























S-1220 


SUBROUTINE AREA (Cont.) 





PN 

-P 1 


PAY(l) 

pi-1 
= ? 

- 1 

PAY<2) 

' 

y w‘ 

- 1 

PAY(3) 


yl 


PAY(4)= 


x PAY (2) 

PAY (3> 
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SUBROUTINE AREA (Cont.) 
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SUBROUTINE BEGIN 
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SUBROUTINE BULK 
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SUBROUTINE COLD 


From MAIN 


Calculate CSTAR 
the PE (static) at 
Nozzle Exit based 
on y= 1.4, then 
thrust(F) and Isp 





Call Subroutine BL 
to obtain Thrust Loss 
Due to Boundary Layer 
(DF) 


PRINT RESULTS 


Subtract Loss 
DP, DVL from 
Theoretical F & 
Isp 


Call Subroutine 
DIV to Find 
Percent Isp Loss 
due to Divergence 
(DVL) 


IS PER 


Calculate 
Additional 
F and Change 
in ISP due to 
Mass Addition 


RETURN 


PRINT 

OUTPUT 
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SUBROUTINE COMB 
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SUBROUTINE CONSRV 



DLS 3 Ajp 
SCH = Pr/Lc 
EX1 = Ay/Ax 



Solve Conservation of 
Momentum equation for 

“n+l,I ( Rb ' <«) 


EXU = b 
EX 12 = b 


n, I * 1/2 
n,I - 1/2 


r ~. r 


EXl = 


Alt 2 
Ax 


YES 



Is 

the flov.- 
field 
axisymmctric / 



EX5 = b /S 

Q, i C 

EX6 = b n , 1+1/2 


EX8 = b n j /Pr 
EX9 = b 0 " 1+1/2 /Pr 

EX7 = b, 1-1/2 /S c 


EX10 = b n>I _ 1/2 /Pr 


EX13 « t!>, Aii 2 b a. 1+1/2. b n, 1-1/2 


EX 14 = & jAii) 2 b n, 1+1/2 _ b n, 1-1/2 

S c S c 


Ax Pr Pr 


EX25 = EX22 - EX24 + EX23 


EX20 = Vl + l/2 

( l -k) 

EX 21 = b n> | - i/j 

k-k) 

P f 2 ~ 

EX 22 = b Q,I + 1/2 


u 2 

“ 23 ■ b 0,I - 1/2 

(1 - l/Pr) — 

2 




EX3o=b nI+1/2 (^.y 

EX31 = \ H/J - s) 


EX36 = ~ EX30 
EX37 = j (EX30 


1+ 1 


EX 32 


EX33 


EX38 = - EX31 
EX39 = EX36 - EX37 + EX38 


] 

EX31[hj+EX 34] 
jhj.j + EX35 J 


(DO Loop) 


EX32 

= ShJofJ 
j 1 14 

l 

EX33 

=Th{ 

tt T 


j 1+1 

I 

EX34 

-fhj , 

«J 


j i-i 

I 

EX35 


“X 
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SUBROUTINE CONSRV (Cont.) 
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SUBROUTINE CONSRV (Cont.) 
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RETURN 


DO 

LOOP 


CONTINUE 


YES 
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SUBROUTINE DYNAM (Cont. ) 


Area of orifice J 
is ADD times input 
cross sectional area, 
CARE A 


Calculate pressure 
ratio and determine 
driving volume , 
NDX3 , and receiving 
volume, NDX4 


Set flow and 
admittance 
equal 0 


fPRES(NDX3)- 
s PRES(NDX4)); 
\ < .001 X 


Calculate unchoked 
admittance, ADMIT, 
and Mass Flow Rate , 
WDOT 


Call FMOLWT 
& SPHEAT to 
determine hot gas 
mol. wt. and y 


NDX3 gas 
v ignited / 


Calculate choked 
admittance, ADMIT, 
and Mass Flow Rate , 
WDOT 


orifice 

choked 



Calculate cold gas 
molecular wt. and 
7, and gas constant 
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SUBROUTINE DYNAM (Cont.) 



Estimate 
(4FL/D), 
Set K = 1 


'CA LL' 
CPRF 
Is duct 
choked 
v ? / 


Calculate 
Reynolds No. , 
REN 


Calculate unchoked 
admittance, mass 
flow (CALL FLOW) 


Calculate choked 
admittance, mass 
flow (CALL CFLOW) 


/REN \ 
/ <2100 \ 
(laminar flow)} 


Calculate laminar 
friction factor 


Use New 
Friction Factor 
K = K + 1 



Determine turbulent 
friction factor 
(CALL FFCF) 


K = 3 


r > 

375 

V J 
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SUBROUTINE DYNA/V1 (Cont. ) 


yis the\. 
f Injection \ 
Model Being 
\ Used ?/ 


1003 


/ Is \ 
NDX4 the 
Combustor. 


Set flow rate from 
the combustor and 
accumulated flow 


/NDX3 one\ 
of the main 
jnanifoldsy 



I Set main flow 
to combustor, WT, 
and main mixture, FM 


NDX3 the 
Combustor 
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SUBROUTINE DYNAM CCont. ) 


Set heat flow rate , 
mass flow rate, and 
sum of admittances 
for NDX3 & NDX4 


/ ' - > 

number oi 

connectors 


Set new mass 
content, new mixture, 
new heat content 
for Volume I 


Loop thru 
volumes 
I = 1 


r T ^ 

Is 

volume 

Ignited 


Determine new 
temperature (CALL TGR) 
and new mol wt 
(CALL FMOLWT) 


Calculate new 
temperature 
and mol wt 


Calculate volume 
capacity and new 
pressure 


Check if sum of 
admittance of Vol. I 
is minimum sum of 
admittance for 
system 


I = 

NNODE 
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DELT “ 
8 CHECK 


f \ 

600 

V J 


DYNAM (Cont. ) 
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SUBROUTINE D.YNAM (Cont. ) 


Save new pressures , 
temperatures , mixtures, 
mass contents , and heat 
flows 


Calculate thrust, 
accumulated impulse , 
integrated Isp 


Print results if 
this is a spark 
arrival time or if 
printout frequency 
dictates 


there 

a 

pilot 
v ? x 


/ the chamber' 
^already been, 
\ lit? / 


tnis a 
spark 
arrival 
v time , 
X ?■ / 


/ * x 

/ pilot lit \ 
but chamber 
V not / 


Set plug 
pressure 
& mixture 
to pilot 
conditions 


Write ignition 
failure message 


/Is\ 

/ the re X. 
/ sufficient X 
potential for 
*\a spark to/ 
\occur/ 


Write spark 
failure 
message 


/ Is the X. 
'spark energyX 
^sufficient for/ 
\ignition/ 
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SUBROL 




EQUILG (LAA, LBB) 



L55 
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SUBROUTINE FG 



Call "FGX" to evaluate the linearized elementary 
reaction equations as a function of temperature, 
pressure, reaction rates, and the existing chemical 
composition, and to evaluate the components of the 
matrix which will be solved to determine the change 
in chemical composition. 


I 


RETURN 
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SUBROUTINE FLUX 
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SUBROUTINE FLUX (Gont.) 
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SUBROUTINE GRAD, NOZ 


From PERF 


For I =1,22 


A (I) = A Log (AR (I)) 



For J = 1,22 
I = 1,6 (corresponding 
to 0F= 4, 6,8 at 100,500 
or 500 and 1000 psia) 
Calculate DMIN (I, J) 


RETURN 


For I = 1, 6 J - 1,22 
DMIN (I,J) = 


which is 


is / dH \ 

UqI 


from NASA 
CR 72601 


From PERF 


Determine Design 
Area Ratio of 
Nozzle 


Generate Table of 
AR(I) from 1=1,22 
Values Range from 
1.001 to 5.0 


RETURN 





1 = 1 , 22 
Calculate 

DADX(I) = f (AR(R) 


L6 


T 




S-1220 


SUBROUTINE GRID 


MINIT = MINIT + 1 
MHALF = 2 x MINIT - 1 



For > A0 j test on difference in properties at lowest 
two grid points to see if a new point should be added 
beneath the grid. Properties examined are velocity, U, 
total enthalpy, H, and specie mass fractions, YSPEC 
(or element mass fractions, ALPHA, for equilibrium 
chemistry option). Any one failure of the test adds the 
grid point. The generalized test is: 

^ NO 

/ Is 


W, 





w 



(DO Loop) 

Transfer information for TELAP(I), H(I), U(I), T(I), 
and ALPHA (K, I), or YSPEC(J,I), HHLD(I),IDA(I), 
THLD(I), and FIX (I), one grid point upward in 
order to add a new point beneath the existing grid. 


MPSI = MPSI + 1 
NPSI = MPSI - 1 
LR = LR + 1 


PO Loop) 

All 0. = - A ip 



12000 , 
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SUBROUTINE GRID (Cont.) 
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SUBROUTINE HEAT (LLA, LLB) 
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Combustor I 

heat transfer Conduction 

characyerisUcs admittance 
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SUBROUTINE HEATT (Cont. ) 
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SUBROUTINE HEATT (Cont. ) 



Find wall 
temperatures 
at this time step 


^ RETURN 

■ 3T 






S-1220 


SUBROUTINE HONC 
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SUBROUTINE HONC (Cont.) 


(DO Loop) 

s '"(i O’-©’ 

' active spec 
Compute DUM - £ / J 

i a W 

all specie 


For the revised Chemical Composition , 
the enthalpy, H, the quadratic enthalpy 
Temperature curve fits, compute the 
New Temperature, T, by inversion. 


DO Loop 
J = 1 to N 


TIME 

= TIME + DT 

AH = 

RHO 


RW 


AT = 

T 

— - 1 . 
TW 



/ WRITE 
/aww, RHO, T, TT , 

FIX (KPSI), DT, TIME, 
TMEF t 


' HXJ (J) < TESTY N 
or 

SAW(J)< TESTY y 


TESTW = 1 .00464 T HXY(J) 1 

TESTX = 1/TESTW L J 

TEST = AW(J)/HXY(J) 


'TEST > TESTW' 
or 

■ TEST < TESTX „ 


INCREMENT J 


/ Is \ 
the DO 
Loop 
complete 
\ ? / 


/ [S \ 

.FIX (KPSI) 
< TMAX 
\ 1 / 


Double Chemistry step size 
DT - 2 x DT 

FIX (KPSI) = 2 x FDC(KPSI) 
NUT = NUT + 1 


/I\ 

temperature 
changing too . 
\. rapidly / 


/ 18 \ 

the density 

changing 

too rapidly 

\ ? / 


FDC(KPSI) = TMAX 


X IB \ 
DT< TMAX 



L<I»J 
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SUBROUTINE HONC (Cont.) 
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SUBROUTINE MARCH 



(DO Loop) 

Compute the Element Mass 
Fractions, ALPHA (K, I), 
from Specie Mass Fractions 
YSPEC (J,I) 


Call "HEAT" 
to compute Temperature, 
C p , and specie enthalpies I 
for all Grid Points 


KOPT = 2 
Specifies 
Enthalpy 
known and 
T to be found 


Call "EQUILC" 
to compute the specie 
Chemical Composition 
for all Grid Points 


Call "DENSE" 


Call "BULK" 

to Compute Density, 
Moleseular Weight 
and Mass Flux for 

w / more than \ - Es 

to Compute the 

r One Grid / r 

bulk Properties 

all Grid Points 

\ Point / 

across the Duct 




Call "VISC" 
to compute 
the Value of 
the Turbulent 
Viscocity Coefficient I 


(DO Loop) 

Compute the Radial Location 
of each streamline, Y (I) 


(DO Loop) 

1 Compute the Vertical Location 
of each streamline, Y (I) 


NO 



J 

L 







(DO Loop) 

Compute the mixing parameter, 

A (1), for all Grid Points, for 
either A zisym metric or 2-D flow 


350 


Call "STEPi^- - 
to compute the Field 


Call J^V’ALL" 
to determine the Wall ' 

Axial Step Size for 
the next Diffusion Step 

i 

\ 

as a Function of the 
Cooling Model Used 
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SUBROUTINE MARCH (Cont.) 





Call "PREPAR" 

A 

1 

r 

u 

to begin the 
finite rate chemistry 
computation for all 
Grid Points 

Y^g y/ Finite Rafce^^v 
Chemistry 
Being Used / 

(DO Loop) 


k 

Set ail WDT ■ . 0 








I no 
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SUBROUTINE MARCH (Cont.) 


/ Are \ 

/ any of the 
velocities or Mass' 
Fractions, across 
the Grid Negative > 
V (DO Loop) / 


WRITE 

Failure Message, 

RX, DX, IOUT, YSPEC, 
RU, RT, RALPHA 


r Is > 

IFINIS = 0 
v ? y 


IFINIS - 1 

MIN IT = IPRESS 

MHALF - 2*M1NIT-1 



X = RX 

n 

n 

P = RP 


PO Loop) 

Adjust arrays of mass 
Fractions, Enthalpy and 
Velocity at the new axial 
Station. Update the 
Streamline elapsed 
residence times. 


Call "FLUX” 

To record the amount of 
mass added to the flow 
field through the wall over 
the last diffusion step 


/ - \ 

transpiration 

or slot wall 

.cooling models > 

\ in use / 

x i / 


/ 1 “ \ 
there ^ 
only one 
grid point 
v ■? / 


Call "GRID" 
to see If additional grid 
points should be added 
to the computation, or 
A"$ 'doubled 
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This subroutine prepares for the solution of the integrated form of the 
chemical kinetic, linearized, ordinary differential equations: 

2 

y (t + At) - y(t ) = I - A At + A 2 ( At) At B 

° 0 2 12 J 
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SUBROUTINE PERF 
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SUBROUTINE PERF (Cont.) 
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SUBROUTINE PERSHF 


From PERF 


¥ 


Calculate Percent 
Recombination /lSP(I,J)j 
For O/F -4,6, 8 at 100 & 
500 or 500 & 1000 psia 
and at the Low & High 
Area Ratios, 20, 40 or 
40 and 100. (J) 



Interpolate to 
^ find ISPK(I) at 
Specified Area 
Ratio 



Set 

ISPK(I) = 0 
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SUBROUTINE PREPAR 



CTM = CTM + FK(I)*U(I) 
KTM = KTM + 1 


If necessary, modify 
Diffusion step size scale 
Factor, XMPS, from 
KTM and CTM so that not 
more than an average of 
XK(5) Chemistry steps 
will be made per diffusion 
step per grid point in future 


YES 


Set IFI 
to tern 
comput 

MIS = 2 
mate 
ation 

1 

r 


RETURN 


/ WRITE Failure / 
r / Message / 
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SUBROUTINE PRESS (PARR, DPARR, X, CGPR, XP, LYY) 




RETURN 


L8 



SUBROUTINE PICRT 


In this subroutine, the solution of 

chemical kinetics matrix, y. = A.. 

i i] 



is obtained: 
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SUBROUTINE PRINT 


g c = !». SOfilij j 

; 

J 

DO LOOP 

pc • - « c , | 

n 


I - 1 to MPSI 


(DO Loop) 

Compute sensible total 
enthalpy ^ 1111(1)) 

equal to 
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SUBROUTINE PRINT (Cont.) 


(DO Loop) WRITE 
*• *V V T V C P{ 
P 1 .y 1 ,Pti) l , w r i 
for I = 1 to MPSI 


(DO I-oop) WRITE 

J. vw^.p,, 

M r V 1 lor y 

I = 1 to MPSI / 


{_ _WRITE 
T, n, W, p, y, M, Ojj, Cp, 

U x , H, oTf, t, (PW) p t , t t . 

5.. P. / 


' there 
only one 
grid 

\ point y 


Pfpsia), XD. y y 

dx > 


(DO Loop, for 1=1 to MPSI) 

WRITE 

, H 0 H 2 0 H„ OH HO, H,0, , 

I. ttj . ttj, aij 2 . aj 2, a , . Oj 2 , 2 . I 


PO Loop, for I = 1, MPSI) 
WRITE 


I. xf. X° X i H 2°, X i H 2, x" 2 . x^.xf 02 . X? 2 ° 2 .l 


/ “X 

'TranspirationN 
wall cooling 
\being uaedy 


/ 18 X 

regenerative > 
wall cooling 
\being used / 

X i / 


/ PO Loop, for I = 1, MPSI) 

/ WRITE 

, Trace Trace H O . J „ • 

I, aj . Xj , aj . “i . O/Fj.Eaj. 


T_,<ou) . F, St, G, 
~ 1 

R*, & t lb / 

e w / 


WRITE 
q w , T w 


/ «• \ 
Slot Wain 

Cooling 
.being used 


C*, T , T , R* T ,C°, tli 1 , 4, 
aw w w p w w 
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SUBROUTINE SI 
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SUBROUTINE SLOT 


Is 

this the 

first axial step \ Y: 
using this / 
cooling 


ES 



Call 'INVERT" 
to Obtain Adiabatic 
Wall Temperature, 
TAW, from HAW 
and A LAW 


Initialise Parameters: 

NSLCR = 1 

. 9 / . 026 t _ 6 / D* \ 

HTCG=A. (- 5 ^) ft (-^) 

FPZ = .001 
IWTLM =* 10 

TWS(K) = T(LR); K - 1 to NS LOT 


.1 


Compute Adiabatic 
Wall Enthalpy, 
HAW From Bulk 
Conditions 


WW Is Wall 

Molecular 

Weight 


(DO Loop) 

Set ALS(J) Equal to 
Chemical Specie at 
the Wall, and 
ALAW(J) to Bulk 
Valves of Species 
Across the Chamber 


DO Loop 
K = NSLCR 
to 

NS LOT 





! [YES 


NSLCR = K-I 

Ij ^ 

N3LCR= NSLOT 


▼ 

i 


Call'"VASC" to 


Evaluate the 


Laminar Viscosity 


at the Wall, VS. 


from Built Conditions 



Compute UB(K), = FFF(K), 
and SF(K) 

(O 1 , F 1 , and I 1 ) 



YES 
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SUBROUTINE SLOKCont.) 


Compute ANW, EMW and CS 
(A*. m, and C») 


/la \ 
MY > 0 
v ? / 


WRITE 

lOUT, NP. r\VTLM, ANW, EMW, 
CS, RRHO, RCP, RVS, TWS (1) 


CALL "VASC" 
to compute coolant 
laminar viscosity, 
VC 


CALL "HOOT" 
to compute coolant 
enthalpy, HCS(K), 
from coolant 
temperature and specie 


DO Loop 
K = 1 to NSLCR 


IW « 0 I — ►{ 140 V-W IW = IW + If 


CALL "VASC" 
to compute laminar 
viscosity, VWS(K), at 
wall from TWS(K) 


Compute RWS(K) and SG(K) 
(P w k and a k ) 


Compute HTC(K), ST(K), 
GGG(K), RYNC(K), and ETA(K) 

(K k , St k , G k . R c k , and T7 k ) 


/ WRITE / / 

/K, IW, KM, HCS(K), VC, / YES / > „ 

VWS(K), TCS(K), TWS(K),/^ \ ? 

RWS(K), SG(K) / \ / 


/ Is \ 
MY > 0 


K, IW, KM, HTC(K), ST(K), GGG(K) 


RHC, RYNC(K), ETALG, ETA(K) 



„ YES, 

/ Is\ 






no] 


YES 


r— TWB 

= TCS(1) <— — < 

iw s i 



? / 



NO | 
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SUBROUTINE SLOT (Cont.) 




EPW _ TWS_(Kj _ 
EPW TWB 1 



WRITE 
IOUT, K, IW, EPW, 
TWB, TWS (K), GRPZ 


TW = TWS (NSLCR) 
MUSE = NUSE + 1 • 


0- 

Increment K 


1 NO 


/ u\ 

YES 

/ the loop 

“ \ 

v on K / 


\complete / 



WRITE 

PROBLEM MESSAGE 
IOUT, ISO BAT, K, TWS(K), 
TW’B, ETA(K), RYNC(K), 
RHC, HTC(K), SG(K) 




Store Information for next atop: 


\ 

f 

RRUT = RUTBL 



RUE = UBL 




RRHO => RHOBL 
RCP = CPBL 
RVS = VS 


1 




RALAW(J) « ALAW{J) 

f YES 
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SUBROUTINE STEP 


/ 18 \ 

'there more 

than one 
grid 

\T/ 


/ - \ 

Finite Uate^ 

Chemistry 
being used 
\ ? / 


Evaluate the step 
size, DX, from 
the velocity and 
the chemistry time 
Step 


f \ 

this the 
first flow 
field step 

v ? y 


Evaluate DX as the smaller 
of either the stroamtube 
radius, or a function of the 
total flowfield length 


Evaluate DX as a function 
of the changes in velocity 
and mass flux 


DLSSCH = Sc(aC) 


/ 13 \ 

/ the flow \ 

field N 

. axisymmetric a 


Evaluate XD as a 
function of center- 
line values for 


Sc, &t£>. p 



(DO Loop) 

Evaluate XD as the minimum 
of all grid point values for 
XD as a function of 
Sc, p , 0, a 


Evaluade XD as 
a function of 
centerline values 
for Sc, I), p 


(DO Loop) 

Evaluate XD as the minimum 
of all grid point values for 
XD as a function of 
Sc, A$ f p ,0, a 


RETURN 


301 mJ DX ' XMPS 


PO Loop) 

Evaluate XD as the minimum 
of XD and the radial distance 
between adjacent grid points 


Set IFINIS = 2 
To terminate 
Computation 


WRITE 

DX, X, XD, XMPS 
















S-1220 


SUBROUTINE TRANSP 



> 

r 

NO 

r 


NED 

* K 





Sct the "edge" conditions for H 
a, it, p, C p , h, andT, uqual 
to the properties of the NED 1 * 1 
grid point. 


aw, 
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S-1220 


SUBROUTINE TRANSP (Cont.) 


0 -, 


Set the "edge" conditions for If , 
• aw 

or,u, p , Cp, h and T, equal to the 
chamber bulk properties. 


■0 




IIIIE 








WRITE 
' NED, ITWLM, NCL.RRUT, 
RUE.RRIIO, RCP.RI1EE 
RHAW, RAWNCL 


T 



Compute FF, WW, WE, WED 

— ► 

Wll = We/Ww 

(F, W , W , W ) 
w e e 


RCW = 1 - o£, 



EKD 


T(LR) 

= TW 





1 



rrw =■ o 

TWB = T(LR) 
TSTAR = TAW 


Cal! "INVERT" 
to compute the temperature | 
TSTAR (T*) from h*. 
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S-1220 


SUBROUTINE 


TRANS P (Cont.) 
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S-1220 


SUBROUTINE ^ 



SC 


Compute the viscosity, XMUT, 
from the HIRSCH Rocket Chamber 
Model 


XMUT is set equal to 
an inputted constant 



S-1220 


SUBROUTINE WALL 


/ this the \ 
~ first axial 
step in the 
.computation 


NUSE -= O 1/;; 
PUT * {£r) 


/“ wallS. VKS 

above the DEER - t 

computational 

\ grid / 


/ is \ 
DEER < .0 


/ WRITE 
PSIW,T(NPSI) 


NPSI -• MFSI - 1 

^ - 

psiw -= psifMPsn 

NO / Is \ 

-A... /l 4 b' - DEERi 

— * •- 

NEW = 1 

PC NT PRNT 


t 

MPSI = MPSI*! 



PSIW = PSl(NPSI) 


LR = MPSI 
JR = NPSI -1 
NED = 1 


, ,,NP , 
' V= (y) v (P~ ), 


Compute adiabatic wall 
total enthalpy, HAW, 
from chamber bulk 
conditions 


Call "L\"VERT" to 
obtain adiabatic wall 
tempersutre - TAW 


RRUSE = RUSE 

. yes /nuse\._ 

SCB = PRB/Le 

NO/ NUSE 

T 

VB - u. 

V 


*7 \ v y*i 



Impose Wall B.C. on 
momentum equation: 


/§Ji.\ = C_ x RRUSE x GAMW 
\h >Ww 2 


/ Is the \ 
T’ranspiratio'Sx YES 

Wall Cooling 
.Model Being / 

Used V 


MPSI 10 >*w 


f \ .^NO x Transpiration^ 
( 300 Wall Cooling 

Vy \ Model in y/ 

(Slot Cooling) \yUse ? / 


' Regeneimtion > 
Cooling 
Impermeable 
\ Model / 
\bsed y 


/ «- \ 

r Issenergetic > 

Impermeable 

s\Vall Model / 

XUsed?/ 


(DO Loop) 

Set element and specie 
.Mass fraction at MPSI 
to values at JR 
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S-1220 


SUBROUTINE WALL (Cont.) 



YES 

Call "EQUILC*’ to compute 
the Chemical Species from 
the Modified Chemical 
Element Composition at the 
Wall and B.C. Grid Points 



Call "HEAT" to 
Compute the Enthaipy 
from the temperature 
and Chemical Specie 
Composition at the Wall 
and B.C. Grid Points 

> P 

2 










FAILURE 



Call "SLOT" to 
Compute the 
Conditions at the 
Slot Cooled Wall 


H(MPSI) = H(JR) 
RH(MPSl) = H{JR) 
KOPT * 2 



Call "HEAT" to 
Evaluate the 
Temperature from 
the Enthalpy at the 
B.C. Grid Point 


JL 

IFINIS = 2 to 







Failure 

Message 


RETURN 


Computation 

W 

W 








25 


Call "DENSE" to 
Compute Density, 
Molecular Weight, 
and Mass Flux at 
the Walt and B.C. 

Grid Points 


Store Information 
for Next Stop: 
RRUSE = RUSE 
RPRB = PRB 
RVB = VB 


RETURN 


©-* 


LR « MPS1 



(Wall is still above 

— » 

RETURN 

computation grid) 
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